- 必要なパッケージのインポート
- テスト用画像生成
- 輪郭抽出
- 中心座標を揃える
- サイズを揃える
- 向きを揃える
- 得られた回転角情報をもとに輪郭の向きを揃える
- 開始点と点列の向きをそろえる
- 本番の楕円フーリエ解析
- 参考文献
前回の記事で楕円フーリエ解析における標準化は結果の解釈を大きく変える可能性があり、きちんと考える必要があることを紹介しました。 また、その中で輪郭に前処理を行い向きや輪郭の開始点を調整しておくと良いこともわかりました。
しかし、前回はシミュレーションの仮定から輪郭の向きや開始点の補正を逆算的に行なっており、この方法は実践的ではありません。 そこで今回はより実践に近い形でのシミュレーションもしていこうと思います。輪郭の向きや開始点の正解がわからないという状況を再現して輪郭補正をやってみたいと思います。
また、今回重要な仮定として、輪郭が線対称となっていることという条件を設けています。この事実から、その軸を元に向きを補正して楕円フーリエ解析を行うという流れで解析例を実験してみたいと思います。実践の中では線対称でない図形を対象にする場合も多いかと思いますが、そのような場合でも、何が図形の向きを揃えるという意味で特徴的な軸がある場合は参考になるかもしれません。
必要なパッケージのインポート
import numpy as np import pandas as pd import cv2 import matplotlib.pyplot as plt import seaborn as sns from pyefd import elliptic_fourier_descriptors, normalize_efd, reconstruct_contour from sklearn.decomposition import PCA from scipy.spatial import KDTree from scipy.optimize import minimize
テスト用画像生成
今回は最初から輪郭の点列が座標として与えられているという状況の一歩前から実験用のテストデータを作成していこうと思います。輪郭の点列データが得られる前は、対象となる物体と背景を識別したニ値画像を持っている状況が多いと思います。そこで今回は前回と同じ輪郭データに対し、ランダムに回転させ、(回転角 )、ランダムに拡大縮小・平行移動し、ランダムな小さいノイズを加えています。この輪郭データを元に
cv2.fillPoly関数を用いて二値画像を作成しています。
あとで図形の回転角psiを統一するという処理が目標の1つになるので、正解の図形の回転角も変数true_psiに保持しておきます。
np.random.seed(ord("🐱")) images = np.zeros((100, 500, 500), np.uint8) i = 0 true_psi = [] for a in np.linspace(-1, 1, 10): for d in np.exp(np.linspace(np.log(1/2), np.log(2), 10)): # テスト用輪郭の生成 theta = np.linspace(0, 2*np.pi, 100)[:-1] x = np.cos(theta) + a * np.cos(theta*2) y = d * np.sin(theta) # ランダムに回転 psi = np.random.uniform(0, np.pi) rot = np.array([[np.cos(psi), -np.sin(psi)], [np.sin(psi), np.cos(psi)]]) x, y = rot @ np.stack((x, y), axis=0) true_psi.append(psi) # 一旦面積を揃えてから平均100の拡大に関する係数をかける tmp = np.stack([x, y], axis=1).reshape((-1, 1, 2)) area = cv2.contourArea(tmp.astype(np.float32)) s = np.random.normal(loc=100, scale=10) x *= s / np.sqrt(area) y *= s / np.sqrt(area) # ランダムなノイズを加える x += np.random.normal(size=len(x)) y += np.random.normal(size=len(y)) # 余白10px含め画像500x500に入らない場合はサイズダウン max_side_size = max(np.ptp(x), np.ptp(y)) if max_side_size > 480: x *= 480 / max_side_size y *= 480 / max_side_size # 余白10px除いた画像中央480x480の領域にランダムに配置されるようにする # または中心座標が200~300に入るようにする x += np.random.uniform(max(10 - x.min(), 200 - x.mean()), min(490 - x.max(), 300 - x.mean())) y += np.random.uniform(max(10 - y.min(), 200 - y.mean()), min(490 - y.max(), 300 - y.mean())) # 輪郭をもとに画像に塗りつぶし pts = np.stack((x, y), axis=1).reshape(-1, 1, 2).round().astype(np.int32) cv2.fillPoly(images[i], [pts], 255, 8, 0) i += 1 # テスト画像のいくつかを抜粋 fig, axes = plt.subplots(1, 5, figsize=(10, 2)) axes = axes.ravel() for i, img in enumerate(images[::20]): ax = axes[i] ax.imshow(img, cmap="gray") plt.show()
輪郭抽出
輪郭抽出の一般的な方法に従いcv2.findContours関数を用いた輪郭検出を行います。描画用の関数plot_contourも用意しておきます。これはtrue_psiには真の回転角、pred_psiには予測された回転角をとり、それぞれ角度の軸を描画します。
def plot_contour(contours, true_psi, pred_psi=None): fig, axes = plt.subplots(10, 10, sharex=True, sharey=True, figsize=(10, 10)) axes = axes.ravel() # 輪郭描画 for ax, cnt in zip(axes, contours): x, y = cnt[:, 0, :].T ax.scatter(x, y, c=range(len(x)), s=1) # 輪郭中心座標の計算 centers = np.array(contours).mean(axis=(1, 2)) # 真の傾き描画 for ax, psi, center in zip(axes, true_psi, centers): ax.axline(center, slope=np.tan(psi), linestyle="dotted", color="gray") # 予測された傾き描画 if pred_psi is not None: for ax, psi, center in zip(axes, pred_psi, centers): ax.axline(center, slope=np.tan(psi), linestyle="dashed", color="orange") for ax in axes: ax.set_aspect(1) #ax.axis('off') return fig, axes
contours = [] centers = [] for img in images: cnts, _ = cv2.findContours(img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE) assert len(cnts) == 1 #contourSamping関数が使えるのであれば等間隔に座標を再サンプリング cnt = cv2.ximgproc.contourSampling(cnts[0], nbElt=100) contours.append(cnt) plot_contour(contours, true_psi) plt.show()
中心座標を揃える
まずは、最も簡単に楕円フーリエ解析で標準化を用いない際に重要な前処理となる中心座標（重心）を揃えることを行いたいと思います。原理的にはどこに揃えてもいいような気はしますが、わかりやすいので原点に合わせる処理を行います。この処理は具体的には各座標の平均をすべての点の座標の値から引いてあげれば良いです。図では原点を青色、輪郭の中心をオレンジ色で描画しており、それらが重なっていることがわかるかと思います。
contours_centered = [cnt - cnt.mean(axis=(0, 1)) for cnt in contours] fig, axes = plot_contour(contours_centered, true_psi) for ax, cnt in zip(axes, contours_centered): ax.plot(0, 0, '*') ax.plot(cnt[:, 0, 0].mean(), cnt[:, 0, 1].mean(), '*') plt.show()
サイズを揃える
楕円フーリエ解析における標準化では第一調和楕円の長軸が一定になるように標準化がされていますが、今回はその情報を用いると細長い輪郭と扁平な輪郭で差が出てしまうことが前回の記事からも予想される懸念としてあるので、面積が1となるような前処理(=面積の平方根で座標の値を割る)を行いたいと思います。この点に関しても、各ケースで揃えるべき部分は変わるかと思います。例えば今回の場合でも矢じりの長さを一定にした場合の輪郭の特徴の差を見たいというケースでは別の前処理が必要になってくるかと思います。こう考えていくと形を揃えるといったことは言葉では簡単に言えますが、なかなか難しいことだと感じます。
contours_resized = [] for cnt in contours_centered: area = cv2.contourArea(cnt) contours_resized.append(cnt / np.float32(np.sqrt(area))) plot_contour(contours_resized, true_psi) plt.show()
向きを揃える
ここで一番の難関、向きを揃えることを目指していこうと思います。もちろんtrue_psiに格納されている情報は実践の場では未知の情報のため使えません。
今回は輪郭の特性として、「線対称である」ということを利用し、対象軸をまず検出して、その角度をもとに輪郭の軸を揃えていきます。
対称性を測るクラスを定義
輪郭の対称性を見つける上で対称性とは何かを定義してあげる必要があります。ここでは輪郭をある軸について線対称になるように変換した輪郭を作り、それと元の輪郭の差異が小さいとき、その軸は対称軸であるだろうと考えます。ここで2つの輪郭の差異はscipy.spatial.KDTreeを用いて測ります。energyメソッドを回転角psiを引数にして呼ぶとその図形中心を通る、x軸とのなす角がpsiの線で対象な図形との差異を計算します。この値が小さいほど、そのpsiは対称軸の回転角らしいと言える指標になります。
class ReflectionSymmetry(): def __init__(self, contour): contour = np.squeeze(contour) self.contour = contour - contour.mean(axis=0) self.tree = KDTree(contour) def reflect(self, psi): n = np.array([-np.sin(psi), np.cos(psi)]) # 法線ベクトル dist = self.contour @ n return self.contour - 2 * dist[:, None] * n def energy(self, psi): reflect_contour = self.reflect(psi) dists, _ = self.tree.query(reflect_contour, k=1) return dists.mean()
輪郭に対してのPCA
広く使われている輪郭の軸の検出で、主成分分析(PCA)があります。点列のx座標とy座標に対してPCAを適用した際の主成分軸が対称線に重なることが多いためです。これはどんなケースでも利用できる方法ではないので、一旦、今回のケースでどうなるかをみてみます。主成分軸が2つ得られるので、そのうちどちらを採用するかは先ほどのenergyメソッドの値が低い方にします。推定された軸をオレンジ色で描画します。
def pred_psi_by_pca(contour): pca = PCA(n_components=2).fit(contour.reshape((-1, 2))) psi1 = np.arctan2(pca.components_[0, 1], pca.components_[0, 0]) psi2 = np.arctan2(pca.components_[1, 1], pca.components_[1, 0]) rs = ReflectionSymmetry(contour) if rs.energy(psi1) < rs.energy(psi2): psi = psi1 else: psi = psi2 return psi pred_psi_pca = [pred_psi_by_pca(cnt) for cnt in contours_resized] fig, axes = plot_contour(contours_resized, true_psi, pred_psi_pca) plt.show()
概ねうまく回転角を推定できていますが、一部大きく外しているものもあります(真の軸(灰色点線)と予測された軸(オレンジ点線)がズレる)。推定が外れるものは比較的縦長でも横長でもない中間的な形状を持っている傾向があります。残念ながら今回は使えなさそう(一部マニュアルで回転角を設定するといった対応が取れる場合はありかも）なので、別の方法を試します。
2次のモーメントを利用
原理的にはPCAと同じ方法ですが、opencvのmoments関数を使った場合も試してみたいと思います。ここで、回転角はこちらのサイトを参考にしました。こちらも推定された回転角がちょうど90どズレる可能性（横長の矢じりの場合）があるので先ほどのPCAと同様の基準で推定された角度とその+90度の2つから1つを選択します。
def pred_psi_by_moment(contour): m = cv2.moments(contour) mu11, mu20, mu02 = m['mu11'], m['mu20'], m['mu02'] psi1 = 0.5 * np.arctan2(2 * mu11, mu20 - mu02) psi2 = (psi1 + (np.pi/2)) % (np.pi*2) rs = ReflectionSymmetry(contour) if rs.energy(psi1) < rs.energy(psi2): psi = psi1 else: psi = psi2 return psi pred_psi_mmn = [pred_psi_by_moment(cnt) for cnt in contours_resized] fig, axes = plot_contour(contours_resized, true_psi, pred_psi_mmn) plt.show()
原理的にはPCAと同じためか、PCAと同様にいくつかの輪郭で推定がうまくいっていないことがわかりました。
EFAの第一調和楕円の軸を利用
前回の楕円フーリエ解析の正規化の記事でも第一調和楕円の長軸の向きを合わせるといった標準化はうまく機能しないことがわかっています。しかし、短軸であれば上手い推定ができているかもしれません。そこで、長軸、短軸のうちどちらが良いかを先ほどと同様にして選択したときにどうなるかもみてみます。第一調和楕円の長軸の回転角はpyefd.normalize_efdで計算可能なのでそちらを利用します。
def pred_psi_by_efa(contour): contour = np.squeeze(contour) efd = elliptic_fourier_descriptors(contour) _, (_, psi1, _) = normalize_efd(efd, return_transformation=True) psi2 = psi1 + (np.pi/2) # 第一調和楕円の短軸の方 rs = ReflectionSymmetry(contour) if rs.energy(psi1) < rs.energy(psi2): psi = psi1 else: psi = psi2 return psi pred_psi_efa = [pred_psi_by_efa(cnt) for cnt in contours_resized] fig, axes = plot_contour(contours_resized, true_psi, pred_psi_efa) plt.show()
残念ながらこれもいくつかの輪郭で推定がうまくいきませんでした。
総当たり法
最終手段、360度を100等分して、そのうち最も対称線らしいものを対称軸として推定します。解像度は3.6度が限界ということになりますが、ここは計算量との兼ね合いで増やすこともできますし、最終的には最適化アルゴリズム(Nelder-Mead法など)を使っても良いと思います。
def pred_psi_by_bruteforce(contour): rs = ReflectionSymmetry(contour) _, psi = min([(rs.energy(psi), psi) for psi in np.linspace(-np.pi/2, np.pi/2, 100)]) return psi pred_psi_bf = [pred_psi_by_bruteforce(cnt.copy()) for cnt in contours_resized] plot_contour(contours_resized, true_psi, pred_psi_bf) plt.plot()
結局これが一番うまくいっているのでこの方法で得られた推定値をもとに図形を逆に回転させ、向きを揃えていきたいと思います。
得られた回転角情報をもとに輪郭の向きを揃える
軸を揃える
先ほど総当たり法で求めた角度と反対方向に回転させるために、回転行列を座標値にかけることで輪郭を回転させます。
contours_rotated = [] for psi, cnt in zip(pred_psi_bf, contours_resized): cnt = np.squeeze(cnt) M = np.array([ [np.cos(psi), np.sin(psi)], [-np.sin(psi), np.cos(psi)], ]) cnt = (M @ cnt.T).T.reshape(-1, 1, 2) contours_rotated.append(cnt) true_psi_rotated = np.array(true_psi) - np.array(pred_psi_bf) plot_contour(contours_rotated, true_psi_rotated) plt.show()
向きをそろえる
まだ矢じりの尖っている方とそうでない方が逆のものも存在するので、それらについて、y軸の分散の値を参考に逆転させます（y軸の分散は尖っている方が小さいという印象のため）。
for cnt in contours_rotated: left = cnt[cnt[:, 0, 0] < 0, 0, 1] right = cnt[cnt[:, 0, 0] > 0, 0, 1] if np.var(left) < np.var(right): cnt[:, 0, 0] *= -1 fig, axes = plot_contour(contours_rotated, true_psi_rotated) plt.show()
開始点と点列の向きをそろえる
最後に、点列の開始点、向き（時計回りか反時計回りか）がを揃えます。具体的にはx軸が最大の点(=矢じりの先端)を開始点になるように、np.rollで要素をずらし、点列の前半と点列の後半のy軸の平均の大小から反時計回りかどうかを検証し、反時計回りになるようにします（反時計回りの場合は先に矢じりの左半分（xy平面上ではy軸の正の方向）に向かうので、y座標の値は前半の方が大きいだろうという予想のもと）。
contours_shifted = [] for cnt in contours_rotated: new_cnt = np.roll(cnt, shift=-cnt[:, 0, 0].argmax(), axis=0) n = cnt.shape[0] former = new_cnt[:n//2, 0, 1] latter = new_cnt[n//2:, 0, 1] if former.mean() < latter.mean(): new_cnt = np.array(list(reversed(new_cnt))) contours_shifted.append(new_cnt) fig, axes = plot_contour(contours_shifted, np.array(true_psi) - np.array(pred_psi_mmn)) plt.show()
contours_final = contours_shifted
本番の楕円フーリエ解析
今回もktchを参考にした楕円フーリエ解析を行います。(https://doc.ktch.dev/notebooks/outline/elliptic_Fourier_analysis.html)。また、今回はShape on R(http://lbm.ab.a.u-tokyo.ac.jp/~iwata/software/shape_r/)のような可視化も行ってみました。これは楕円フーリエ記述子の値を主成分ごとに分けて一個ずつ変化させていった時にどのような輪郭が得られるかを重ねて描画するものです。
coefs = [elliptic_fourier_descriptors(cnt[:, 0, :]).flatten() for cnt in contours_final] pca = PCA(n_components=5) coefs_transformed = pca.fit_transform(coefs) plt.bar(range(5), pca.explained_variance_ratio_) plt.show()
まずは楕円フーリエ記述子で記述される今回の解析対象の輪郭はおおよそ2~3主成分で説明可能なことがわかります。 それでは、ktchのNotebookにあるようなMorphospaceの描画と、Shape on Rにあるような各主成分が説明する輪郭の変動の描画を行ってみます。
def get_pc_scores_for_morphospace(ax, num = 5): xrange = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1],num) yrange = np.linspace(ax.get_ylim()[0], ax.get_ylim()[1],num) return xrange, yrange # 参考：https://doc.ktch.dev/notebooks/outline/elliptic_Fourier_analysis.html def plot_recon_morphs(pca, fig, ax, n_PCs_xy=[1, 2], morph_num=5, morph_scale=0.8, morph_color="gray", morph_alpha=0.7): pc_scores_h, pc_scores_v = get_pc_scores_for_morphospace(ax, morph_num+2) pc_scores_h, pc_scores_v = pc_scores_h[1:-1], pc_scores_v[1:-1] for pc_score_h in pc_scores_h: for pc_score_v in pc_scores_v: pc_score = np.zeros(pca.n_components_) n_PC_h, n_PC_v = n_PCs_xy pc_score[n_PC_h-1] = pc_score_h pc_score[n_PC_v-1] = pc_score_v arr_coef = pca.inverse_transform([pc_score]) ax_width = ax.get_window_extent().width fig_width = fig.get_window_extent().width fig_height = fig.get_window_extent().height morph_size = morph_scale*ax_width/(fig_width*morph_num) loc = ax.transData.transform((pc_score_h, pc_score_v)) axins = fig.add_axes([loc[0]/fig_width-morph_size/2, loc[1]/fig_height-morph_size/2, morph_size, morph_size], anchor='C') coords = reconstruct_contour(arr_coef[0].reshape(-1, 4)) x = coords[:,0] y = coords[:,1] axins.plot(x.astype(float),y.astype(float),color=morph_color, alpha = morph_alpha) axins.axis('equal') axins.axis('off') # 参考：http://lbm.ab.a.u-tokyo.ac.jp/~iwata/software/shape_r/ def plot_pca_recon(pca, i, unit, normalize=False, ax=None): pcscores = np.zeros((3, pca.n_components_)) pcscores[:, i] = np.array([-1, 0, 1]) * unit * np.sqrt(pca.explained_variance_[i]) coeffs = pca.inverse_transform(pcscores) coeffs = coeffs.reshape((3, -1, 4)) contours = [reconstruct_contour(coeff) for coeff in coeffs] ax.plot(contours[0][:, 0], contours[0][:, 1], 'C0--') ax.plot(contours[1][:, 0], contours[1][:, 1], 'C2') ax.plot(contours[2][:, 0], contours[2][:, 1], 'C3--') ax.set( aspect=1, title=f"PC{i+1}" ) return ax
import pandas as pd df_pca = pd.DataFrame(coefs_transformed, columns=[f"PC{i+1}" for i in range(5)]) fig, ax = plt.subplots() sns.scatterplot(data=df_pca, x="PC1", y="PC2", ax = ax, legend = True) plot_recon_morphs(pca, fig=fig, ax=ax) plt.show() fig, axes = plt.subplots(1, 3, figsize=(9, 3)) for i in range(3): plot_pca_recon(pca, i, 1, normalize=False, ax=axes[i])
前回の記事と同様にいい感じで矢じりの扁平/細長と尻部分の窪みの大小が主成分として抽出できたのではないでしょうか。これで、楕円フーリエ解析で関数備え付けの標準化(Normalization)が機能しない場合でもこのような前処理を行なっていくことで解析が可能であるという道筋が見えてきました。向きの統一はかなり難しい点で対応方法もケースバイケースだと思うので、最近ではAIの力を使うなりしたいところだと思いました(実際今回の記事を書く上でも生成AIに相談していたりします)。
