楕円フーリエ解析と標準化について

楕円フーリエ解析(Elliptic Fourier Analysis, EFA) は閉曲線の形状解析で広く用いられている手法です1。この解析で得られる係数は、輪郭の回転やスケールなどに対して不変ではありません。例えば形が画像中で回転しても解析者は同じ形とみなすことが多いため、異なる向きの同じ輪郭を入力したときに係数に差異が生じるのは不都合です。そのため、多くのライブラリでは標準化(Normalization)*1のオプションが提供されており、第1調和楕円を基準にして輪郭の向きやサイズを揃える処理が行われます1,2,3。

しかし一方で、標準化によって意図していた結果が得られなくなってしまうといった可能性も指摘されています4。そこで今回は単純なケースについてシミュレーションで実験してその状況についてと解決のアイデアについてみていきたいと思います。

必要なパッケージのロード

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from itertools import product
from pyefd import elliptic_fourier_descriptors, reconstruct_contour
from ktch.outline import EllipticFourierAnalysis
from sklearn.decomposition import PCA

テスト用の輪郭生成

今回は2つのパラメータ と を変化させ、矢じりのような形状の集合を生成してみました。以下、輪郭データの生成と描画のコードになります。各輪郭は100個の点の座標列で、重心は原点に揃えています。

# 輪郭の差異を出すためのパラメータ
a2 = np.linspace(-1, 1, 10)
d1 = np.exp(np.linspace(np.log(1/2), np.log(2), 10))

# 輪郭データの生成
contours = []
for a2_i, d1_i in product(a2, d1):
    # 弧長を計算して等間隔になるようにthetaの値を決める(下準備)
    theta = np.linspace(0, 2*np.pi, 1000)
    dtheta = np.diff(theta)
    theta = theta[:-1]
    dx = -2 * a2_i * np.sin(theta*2) - np.sin(theta)
    dy = d1_i * np.cos(theta)
    dt = np.sqrt(dx**2 + dy**2) * dtheta
    t = np.cumsum(dt)
    t = np.concatenate(([0], t))
    target = np.linspace(0, t[-1], 101)
    idx = np.searchsorted(t, target, side='left')[:-1]
    theta = theta[idx]
    # 実際の輪郭生成はここから
    x = np.cos(theta) + a2_i * np.cos(theta*2)
    y = d1_i * np.sin(theta)
    contour = np.stack((x, y), axis=1)
    contour -= contour.mean(axis=0) # 輪郭重心の位置を(0, 0)にする
    contours.append(contour)

# 全100データを10x10で表示
fig, axes = plt.subplots(10, 10)
for i, (contour) in enumerate(contours):
    ax = axes[i//10, i%10]
    ax.scatter(contour[:, 0], contour[:, 1], c=np.linspace(0, 1, 100), cmap='viridis', s=.5)
    ax.set_aspect(1)
    ax.axis('off')
plt.show()

テスト用輪郭

上記コードにより、色で点列の進行方向が示された矢じり型の輪郭が描画されました（紺色が開始点、点列が進むほど黄色に近づきます）。以降の解析ではこれら各100個の点で表現された100個の輪郭を対象に楕円フーリエ解析を行い、得られた係数についての主成分分析(PCA)の結果を可視化していきます。

楕円フーリエ解析 with ktch

今回用いたのは"ktch"と呼ばれるPythonパッケージです。

このパッケージのExample notebooksにわかりやすく解析の流れが例示されていてそこで定義されている関数を参考に楕円フーリエ記述子についての主成分分析で得られる第1～第3主成分までの散布図とその散布図上の均等にサンプリングした点から再構築した輪郭、さらに、第5主成分までの寄与率の図を一気に描画する補助関数を定義しておきます。

def plot_pca(coef=None, pca=None, efa=None):
    pcscore = pca.fit_transform(coef)

    # pcaの散布図と寄与率の棒グラフ
    fig, axes = plt.subplots(2, 2, figsize=(10, 8))

    for i in range(2):
        for j in range(2):
            ax = axes[i, j]
            if j > i:
                ax.bar([f"PC{i+1}" for i in range(5)], pca.explained_variance_ratio_[:5])
                continue
            ax.scatter(pcscore[:, j], pcscore[:, i+1], s=8)
            ax.set_xlabel(f"PC{j+1}")
            ax.set_ylabel(f"PC{i+2}")
            

            # 主成分が変動したときの輪郭の再構築
            xmin = pcscore[:, j].min()
            xmax = pcscore[:, j].max()
            ymin = pcscore[:, i+1].min()
            ymax = pcscore[:, i+1].max()
            xs = np.linspace(xmin, xmax, 5)
            ys = np.linspace(ymin, ymax, 5)
            ax.set_xlim(xmin - (xmax - xmin) * 0.2, xmax + (xmax - xmin) * 0.2)
            ax.set_ylim(ymin - (ymax - ymin) * 0.2, ymax + (ymax - ymin) * 0.2)
            for x in xs:
                for y in ys:
                    pcscore_manual = np.zeros((1, pca.n_components_))
                    pcscore_manual[:, j] = x
                    pcscore_manual[:, i+1] = y
                    coef_manual = pca.inverse_transform(pcscore_manual)
                    ax_width = ax.get_window_extent().width
                    fig_width = fig.get_window_extent().width
                    fig_height = fig.get_window_extent().height
                    morph_size = .8 * ax_width / (fig_width * 5)
                    loc = ax.transData.transform((x, y))
                    axins = fig.add_axes([loc[0]/fig_width-morph_size/2, loc[1]/fig_height-morph_size/2, morph_size, morph_size], anchor='C')
                    contour_manual = efa.inverse_transform(coef_manual)[0]
                    axins.plot(contour_manual[:, 0], contour_manual[:, 1], color="black", alpha=0.8, linewidth=.8)
                    axins.axis('equal')
                    axins.axis('off')           
    plt.plot()

標準化あり(norm=True)の結果

まずは、ktch.outline.EllipticFourierAnalysisを用いて標準化あり(norm=True)の楕円フーリエ解析を行い、得られた係数行列に主成分分析を適用してその結果を先ほどの補助関数を使ってプロットしてみます。

efa = EllipticFourierAnalysis(n_harmonics=40)
pca = PCA(n_components=None)
coef = efa.fit_transform(np.array(contours), norm=True)

plot_pca(coef, pca, efa)

ktch標準化あり

一見うまくいっているように見えますが、PC1とPC2の散布図をよく見るとPC1軸が矢じりの左右方向、PC2軸が上下方向を区別していることがわかります。しかし、形の解析において、矢じりの向きは、輪郭を画像から抽出する際にたまたま右向きに写っていただけでそれらは同じ形として認識させたい場合が多いかと思います。

また、そもそも入力データには矢じりの向きは左右しかなかったのになぜ上下方向も区別するような結果となっているのでしょうか？この点の検証のため、「標準化あり」で計算された楕円フーリエ記述子からの輪郭の再構築を行ってみたいと思います。

fig, axes = plt.subplots(10, 10)
re_contours = efa.inverse_transform(coef)
for i, (contour, re_c) in enumerate(zip(contours, re_contours)):
    ax = axes[i//10, i%10]
    ax.plot(contour[:, 0], contour[:, 1])
    ax.plot(re_c[:, 0], re_c[:, 1])
    ax.set_aspect(1)
    ax.axis('off')
plt.show()

標準化ありの再構築輪郭

すると、標準化によって輪郭の向きが変更されていることがわかります(入力輪郭:青線、再構築輪郭:オレンジ線)。楕円フーリエ解析における一般的な標準化のプロセスの一つとして、先ほども述べた通り、第1調和楕円の向きをそろえるように標準化するといったことがされます。その結果、縦長の輪郭は横向きに倒されていることがわかります（上図６列目の一部と7列目以降）。

ということで今回のケースでは関数備え付けの標準化(norm=True)では余計な処理が入ってしまうことがわかったので、Falseにして再解析してみます。

標準化なし(norm=False)の結果

標準化をオフにして、同じEFAおよびPCAを実行します。

efa = EllipticFourierAnalysis(n_harmonics=40)
pca = PCA(n_components=None)
coef = efa.fit_transform(np.array(contours), norm=False)

plot_pca(coef, pca, efa)

ktch標準化なし

無事、矢じりの上下向きの違いを代表するような主成分がなくなりました（一回目の解析のPC2）。実際、得られた記述子を用いて輪郭の再構築を行うと以下のようになります。

fig, axes = plt.subplots(10, 10)
re_contours = efa.inverse_transform(coef)
for i, (contour, re_c) in enumerate(zip(contours, re_contours)):
    ax = axes[i//10, i%10]
    ax.plot(contour[:, 0], contour[:, 1])
    ax.plot(re_c[:, 0], re_c[:, 1])
    ax.set_aspect(1)
    ax.axis('off')
plt.show()

標準化なし輪郭再構築

入力輪郭(青線)と再構築輪郭(オレンジ線)がほぼ重なっており、標準化の影響がないことがわかります。ところで、今回の入力輪郭には左右を向いた矢じりが混在しているため、norm=FalseでもPC1が向き(右向きvs左向き)の軸になっています。この向きの差を、形状の差として解析上区別したい場合はこのままでもよいのですが、向きに意味を持たない場合は望ましくありません。

輪郭の向き補正+標準化なしの結果

そこで、輪郭の向きをあらかじめ統一してみようと思います。具体的には最初の50個の輪郭を左右反転し(x座標に-1を掛け)、すべて右向きにそろえます。実解析では輪郭の向きを数学的に推定して揃えるのは難しいですが、今回はシミュレーションデータなので簡単に処理できます。補正後の輪郭をプロットすると、すべて同じ方向を向いた矢じりになっていることが確認できます。

# 向きをそろえた輪郭データ作成
contours2 = []
for i, contour in enumerate(contours):
    sign = -1 if i < 50 else 1
    x = sign * contour[:, 0]
    y = contour[:, 1]
    contours2.append(np.stack((x, y), axis=1))

fig, axes = plt.subplots(10, 10)
for i, (contour) in enumerate(contours2):
    ax = axes[i//10, i%10]
    ax.scatter(contour[:, 0], contour[:, 1], c=np.linspace(0, 1, 100), cmap='viridis', s=.5)
    ax.set_aspect(1)
    ax.axis('off')
plt.show()

向き補正後の輪郭

無事矢じりの向きがそろったところで、再びEFA+PCA解析を行ってみます。

efa = EllipticFourierAnalysis(n_harmonics=40)
pca = PCA(n_components=None)
coef = efa.fit_transform(np.array(contours2), norm=False)

plot_pca(coef, pca, efa)

ktch、向き補正+標準化なし

今度は確かにPC1が矢じりの向きを表す成分ではなくなり、矢じりの向きによるばらつきが消えました。しかし、散布図をよく見ると、再構築される輪郭にPC1軸上で対称性が見られます。すなわち、主成分の値の符号が逆になるだけで形状が変わらない領域が残っており、何かまだ不適切な扱いがあることが示唆されました。

輪郭の向き補正+開始点補正+標準化なしの結果

楕円フーリエ解析の結果に影響を与える別の要因として、輪郭の開始点があります1。本例では輪郭を色付きでプロットしており(先頭が紺→末尾が黄色)、輪郭の上半分(左右反転した50個)を見ると、矢じりの頂点側ではなく尻尾側を始点にしていることが分かります。言い換えれば、左右反転群の開始点が矢じりの先端側でない状態でした。今回は簡易的に輪郭長の半分（=50点）だけ点列をロールシフトすることで開始点をそろえてみます。

contours3 = []
for i, contour in enumerate(contours):
    if i < 50:
        x = -np.roll(contour[:, 0], shift=50) # 100の点で表現されているので50ずらすとスタート地点が半分ずれる
        y = -np.roll(contour[:, 1], shift=50)
    else:
        x = contour[:, 0]
        y = contour[:, 1]
    contours3.append(np.stack((x, y), axis=1))

fig, axes = plt.subplots(10, 10)
for i, (contour) in enumerate(contours3):
    ax = axes[i//10, i%10]
    ax.scatter(contour[:, 0], contour[:, 1], c=np.linspace(0, 1, 100), cmap='viridis', s=.5)
    ax.set_aspect(1)
    ax.axis('off')
plt.show()

向き、開始点補正後の輪郭

これで「向き」に加えて、「開始点」(色で表現)も統一されました。改めてEFAとPCAを実行します。

efa = EllipticFourierAnalysis(n_harmonics=40)
pca = PCA(n_components=None)
coef = efa.fit_transform(np.array(contours3), norm=False)

plot_pca(coef, pca, efa)

ktch、向き・開始点補正後、標準化なし

この結果では、PC1とPC2が形状の本質的な変動成分を表すようになっています。たとえばPC1の値が大きい方へ移動させると矢じりがより扁平になり、小さい方では尖った細長い矢じりになります。またPC2では、大きい側がカッコいい矢じりっぽい形、小さい側はかわいらしいおにぎり型に近づく様子が観察されました。これらの主成分軸は、最初に設定したパラメータとの影響と概ね対応していることが分かります。なおktchでは直流成分(DC成分)として,も結果の返り値に含まれるため、輪郭の重心位置が異なるデータでは慣習的に標準化が必要ですが、今回は全て重心を原点にそろえているため特に問題は生じていません。標準化に関しては他にも留意点があり、最近の研究でも議論されています6。

pyefdパッケージによる検証

同様の解析をPythonの別ライブラリpyefdでも行ってみました。pyefdもKuhl & Giardina(1982)の手法を実装しており、normalize=Trueオプションで同様の標準化を行います7。以下の関数はPCA結果を図示する補助関数で、pyefd.elliptic_fourier_descriptorsの出力を利用します。

def plot_pca_pyefd(coef=None, pca=None):
    pcscore = pca.fit_transform(coef)

    # pcaの散布図と寄与率の棒グラフ
    fig, axes = plt.subplots(2, 2, figsize=(10, 8))

    for i in range(2):
        for j in range(2):
            ax = axes[i, j]
            if j > i:
                ax.bar([f"PC{i+1}" for i in range(5)], pca.explained_variance_ratio_[:5])
                continue
            ax.scatter(pcscore[:, j], pcscore[:, i+1], s=8)
            ax.set_xlabel(f"PC{j+1}")
            ax.set_ylabel(f"PC{i+2}")
            ax.set_xlim(-10, 10)

            # 主成分が変動したときの輪郭の再構築
            xmin = pcscore[:, j].min()
            xmax = pcscore[:, j].max()
            ymin = pcscore[:, i+1].min()
            ymax = pcscore[:, i+1].max()
            xs = np.linspace(xmin, xmax, 5)
            ys = np.linspace(ymin, ymax, 5)
            ax.set_xlim(xmin - (xmax - xmin) * 0.2, xmax + (xmax - xmin) * 0.2)
            ax.set_ylim(ymin - (ymax - ymin) * 0.2, ymax + (ymax - ymin) * 0.2)
            for x in xs:
                for y in ys:
                    pcscore_manual = np.zeros((1, pca.n_components_))
                    pcscore_manual[:, j] = x
                    pcscore_manual[:, i+1] = y
                    coef_manual = pca.inverse_transform(pcscore_manual)
                    ax_width = ax.get_window_extent().width
                    fig_width = fig.get_window_extent().width
                    fig_height = fig.get_window_extent().height
                    morph_size = .8 * ax_width / (fig_width * 5)
                    loc = ax.transData.transform((x, y))
                    axins = fig.add_axes([loc[0]/fig_width-morph_size/2, loc[1]/fig_height-morph_size/2, morph_size, morph_size], anchor='C')
                    contour_manual = reconstruct_contour(coef_manual.reshape(-1, 4))
                    axins.plot(contour_manual[:, 0], contour_manual[:, 1], color="black", alpha=0.8, linewidth=.8)
                    axins.axis('equal')
                    axins.axis('off')           
    plt.plot()

coef = []
pca = PCA(n_components=None)
coef = np.array([
    elliptic_fourier_descriptors(contour, order=40, normalize=True).ravel()
    for contour in contours
])
plot_pca_pyefd(coef, pca) 

pyefd、標準化あり

coef = []
pca = PCA(n_components=None)
coef = np.array([
    elliptic_fourier_descriptors(contour, order=40, normalize=False).ravel()
    for contour in contours
])
plot_pca_pyefd(coef, pca) 

pyefd、標準化なし

coef = []
pca = PCA(n_components=None)
coef = np.array([
    elliptic_fourier_descriptors(contour, order=40, normalize=False).ravel()
    for contour in contours3
])
plot_pca_pyefd(coef, pca) 

pyefd、向き・開始点補正、標準化なし

R/Momocsによる検証

最後に、R言語のMomocsパッケージ8でも同様の解析を行いました。最初にPythonで生成した輪郭をCSVファイルに保存し、Momocsのefourier()で解析しました。以下は、CSV保存のPythonコードと、実際に実行したRスクリプトの抜粋です。

# 最初の輪郭情報の保存
pd.concat([pd.DataFrame(contour, columns=["x", "y"]).assign(i=i) for i, contour in enumerate(contours)]).to_csv("contours.csv")
# 向き、スタート地点補正済みの輪郭情報の保存
pd.concat([pd.DataFrame(contour, columns=["x", "y"]).assign(i=i) for i, contour in enumerate(contours3)]).to_csv("contours3.csv")

library(Momocs)

png("momocs_fig1.png", width = 600, height = 400, pointsize = 20)
read.csv("contours.csv") |>
  dplyr::select(i, x, y) |>
  tidyr::nest(coord = -i) |>
  dplyr::pull(coord) |> 
  purrr::map(as.matrix) |>
  Out() |>
  efourier(40, norm=T) |>
  PCA() |>
  plot_PCA(axes = c(1, 2))
dev.off()

png("momocs_fig2.png", width = 600, height = 400, pointsize = 20)
read.csv("contours.csv") |>
  dplyr::select(i, x, y) |>
  tidyr::nest(coord = -i) |>
  dplyr::pull(coord) |> 
  purrr::map(as.matrix) |>
  Out() |>
  efourier(40, norm=F) |>
  PCA() |>
  plot_PCA(axes = c(1, 2))
dev.off()

png("momocs_fig3.png", width = 600, height = 400, pointsize = 20)
read.csv("contours3.csv") |>
  dplyr::select(i, x, y) |>
  tidyr::nest(coord = -i) |>
  dplyr::pull(coord) |> 
  purrr::map(as.matrix) |>
  Out() |>
  efourier(40, norm=F) |>
  PCA() |>
  plot_PCA(axes = c(1, 2))
dev.off()

• 関数備え付けの標準化

Momocs、標準化あり

• 標準化なし

Momocs、標準化なし

• 輪郭補正後標準化なし

Momocs、向き・開始点補正、標準化なし

まとめ

1. 楕円フーリエ解析では標準化が重要
2. 関数備え付けの標準化ではうまくいかないときがある
3. うまくいかないときは輪郭の前処理で向きや開始点などをそろえて標準化なしで解析する方法がある
4. ktch, pyefd, Momocsで概ね同様の結果が得られる

参考文献

1. Kuhl, Frank P., and Charles R. Giardina. "Elliptic Fourier features of a closed contour." Computer graphics and image processing 18.3 (1982): 236-258.
2. https://www.naro.affrc.go.jp/archive/niaes/sinfo/result/result11/result11_28.html
3. http://lbm.ab.a.u-tokyo.ac.jp/~iwata/software/shape_r/shape_r_ppt.pdf
4. https://momx.github.io/Momocs/reference/efourier.html
5. https://doc.ktch.dev/index.html
6. Wu, Hui, et al. "Reliable and superior elliptic Fourier descriptor normalization and its application software ElliShape with efficient image processing." arXiv preprint arXiv:2412.10795 (2024).
7. https://github.com/hbldh/pyefd
8. https://momx.github.io/Momocs/

*1:個人的にはNormalizationの訳語としては「正規化」の方がしっくりくるのですが、多くの日本語文献で「標準化」として訳されているので倣います。