簡単な輪郭（形）の解析をしてみました

輪郭（形）というと一見定量的に評価するのが難しそうな対象ですが、実は意外とできるということで簡単に実験してみました。

今回の方法の基盤となるのはフーリエ記述子という概念です。

解析対象は最近スーパーでニンジンを大量に買ったのでそれを使いました。

以下に解析手順の全体像を示します。

必要なライブラリのインポート

import numpy as np
import cv2
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
from pathlib import Path
from scipy.cluster.hierarchy import linkage, dendrogram
from sklearn.decomposition import PCA

import torchvision.transforms.functional as F

from torchvision.utils import make_grid
from torchvision.io import read_image
from torchvision.models.detection import maskrcnn_resnet50_fpn_v2, MaskRCNN_ResNet50_FPN_V2_Weights
from torchvision.utils import draw_segmentation_masks

解析対象の画像を読み込んで表示

輪郭をデータとして解析する上では画像から解析するのがおそらく最も便利です。 そこで、今回買ったニンジンをスマホで３パターン撮りました。 （ニンジンはこの後スタッフが美味しくいただきました。）

filenames = [
    str(Path('assets') / 'carrot1.jpeg'),
    str(Path('assets') / 'carrot2.jpeg'),
    str(Path('assets') / 'carrot3.jpeg'),
]
n_files = len(filenames)
images = [read_image(filename) for filename in filenames]

plt.imshow(F.to_pil_image(make_grid(images)))
plt.axis('off')
plt.show()

解析対象のマスク画像を取得

最初にやるべきことは解析対象のマスク画像（ここでは個体を識別しつつ、ニンジンが写っているピクセルで1、背景が写っているピクセルで0をとるような二値画像）の取得です。

今回はMask R-CNN(https://arxiv.org/abs/1703.06870)を用いました。PyTorchには学習済みの重みが存在するのでそちらを使うことにします。このモデルは様々な物体の検出に用いることができ、その中にはニンジンも含まれているので使っています。

結果を描画してみるとうまくニンジンと背景を個体別に識別できていることがわかります。

weights = MaskRCNN_ResNet50_FPN_V2_Weights.DEFAULT
transforms = weights.transforms()

transformed_images = [transforms(image) for image in images]

model = maskrcnn_resnet50_fpn_v2(weights=weights, progress=False)
model = model.eval()

# 予測
outputs = model(transformed_images)
masks = [output['masks'] for output in outputs]

proba_threshold = .7
score_threshold = .85

boolean_masks = [
    output['masks'][output['scores'] > score_threshold] > proba_threshold
    for output in outputs
]

# 描画
img_with_masks = [
    draw_segmentation_masks(image, boolean_mask.squeeze(1))
    for image, boolean_mask in zip(images, boolean_masks)
]
plt.imshow(F.to_pil_image(make_grid(img_with_masks)))
plt.axis('off')
plt.show()

輪郭情報の抽出

マスク画像ができたら輪郭の抽出は単純でOpenCVのfindContours関数をかけるだけです。

ここで輪郭のデータ構造について軽く説明してみます。輪郭は数学的にはxy平面において媒介変数tを用いて

といった形式で表現できます。(例えば円周であれば、f,gがそれぞれcos,sinになります。)

一方で画像解析で扱う輪郭を式で表現するのは難しいので、tについていい感じに標本化して得られる座標の集合として表します。

つまり、輪郭をなぞったときのピクセルの位置を順番に集めたものとして扱われます。

contours_list = []
for boolean_mask in boolean_masks:
    contours = []
    for mask in boolean_mask:
        mask_np = mask.squeeze(0).numpy().astype(np.uint8)
        con, _ = cv2.findContours(mask_np, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
        contours.append(con[0])
    contours_list.append(contours)

# 描画
plt.imshow(cv2.hconcat([
    cv2.drawContours(img.permute(1, 2, 0).numpy().copy(), contours, -1, (0, 255, 255), 30)
    for img, contours in zip(images, contours_list)
]))
plt.axis('off')
plt.show()

輪郭を重ね合わせるための情報を計算する

先ほど示した輪郭の数学的な表現形式を見てみるとtについて周期関数と見ることができ（一周すると元に戻ってくるので）、フーリエ級数展開が適用できます。

実際には離散データなので離散フーリエ変換を行うことになりますがこの時に得られる数列のことをフーリエ記述子と呼び、これが輪郭の情報をより規格的に示していることになります。

これにより、様々に複雑な形を示す輪郭を数値的に評価・比較できるようになるということです。

ここではものすごく端折って説明したので多少不正確なことがあるかもしれません。フーリエ記述子については他にもたくさんの記事が出ているので参考にしてみてください。

2つの輪郭から得られたこの数列の差を最小にするような変換（平行移動・回転・拡大縮小）を計算することで輪郭を重ねることができるようになります。 今回は写真に様々な角度でニンジンが写っているので、写っているニンジンの1つを参照としてそれに重ねるのではなく、ニンジンの形を模した逆三角形に重なるような変換を計算しています。

fit = cv2.ximgproc.createContourFitting()
transformed_contours_list = []
transformed_contours_without_scale_list = []
indv_imgs_list = []

# ニンジンの形を模した逆三角形に輪郭をフィッティングすることで整列された輪郭を得る
carrot_shape_model = np.array([[[0, 0]], [[100, 650]], [[200, 0]]], dtype=np.int32)
for img, contours in zip(images, contours_list):
    img_np = img.permute(1, 2, 0).numpy()
    n_contours = len(contours)

    transformed_contours = []
    transformed_contours_without_scale = []
    indv_imgs = []
    for i in range(n_contours):
        # 重ね合わせに必要な変換の計算
        alphaPhiST, dist = fit.estimateTransformation(contours[i], carrot_shape_model)
        dst = cv2.ximgproc.transformFD(contours[i], alphaPhiST, fdContour=False)
        alphaPhiST2 = alphaPhiST.copy()

        # 拡大縮小しない版も計算しておく
        alphaPhiST2[0, 2] = 1.0
        dst2 = cv2.ximgproc.transformFD(contours[i], alphaPhiST2, fdContour=False)
        transformed_contours.append(dst)
        transformed_contours_without_scale.append(dst2)

        # 個体ごとのサムネイル的なものを作成
        center = contours[i].mean(axis=0).squeeze(0)
        phi_degree = -alphaPhiST[0, 1] * 180 / np.pi
        s = alphaPhiST[0, 2]
        M = cv2.getRotationMatrix2D(center, phi_degree, s)
        M[:, 2] += alphaPhiST[0, 3:] +40
        img_conv = cv2.warpAffine(img_np, M, (280, 750))
        cv2.drawContours(img_conv, [(dst+40).astype(np.int32)], 0, (0, 255, 100), 5)
        indv_imgs.append(img_conv)
        
    transformed_contours_list.append(transformed_contours)
    transformed_contours_without_scale_list.append(transformed_contours_without_scale)
    indv_imgs_list.append(indv_imgs)

輪郭を重ね合わせて描画する

輪郭を重ね合わせるのに必要な変換が計算できたので実際に変換された輪郭を描画してみます。

こうすることで各個体の形の違いの程度といったものが見えてくるかと思います。

def contour_plot(contours, ax):

    colors = [[int(c*255) for c in col] for col in mpl.colormaps['Dark2'].colors]

    # 全ての輪郭含むようなBounding Boxを計算
    xmin, ymin = np.array([cnt.min(axis=(0, 1)) for cnt in contours]).min(axis=0)
    xmax, ymax = np.array([cnt.max(axis=(0, 1)) for cnt in contours]).max(axis=0)
    w = xmax - xmin
    h = ymax - ymin
    
    # Bounding boxの幅が200pxになり、余白20pxを含むようにする
    ratio = 200 / w 
    mergin = 20
    nw, nh = int(w*ratio+20), int(h*ratio+20)
    canvas = np.zeros((nh, nw, 3), dtype=np.uint8)
    canvas.fill(255)

    for j, cnt in enumerate(contours):
        # Bounding box 似合わせスケーリングされたcanvasに合うように輪郭もスケーリングしてから描画
        resized_cnt = (cnt - np.array([xmin, ymin])) * ratio + mergin 
        cv2.drawContours(canvas, [resized_cnt.astype(np.int32)], 0, colors[j%len(colors)], 2)
    ax.imshow(canvas)
    ax.axis('off')

fig, axs = plt.subplots(ncols=n_files, squeeze=False)
for i, contours in enumerate(transformed_contours_list):
    contour_plot(contours, axs[0, i])
    axs[0, i].set_ylim(800, 0)
    axs[0, i].set_title(filenames[i])

plt.show() 

ちなみに変換のうち、平行移動と回転だけ許し、拡大縮小をしないで重ね合わせると以下のようになります。 ニンジンでは横幅よりも長さにより違いがみられるのがわかるかと思います。

fig, axs = plt.subplots(ncols=n_files, squeeze=False)
for i, contours in enumerate(transformed_contours_without_scale_list):
    contour_plot(contours, axs[0, i])
    axs[0, i].set_ylim(800, 0)
    axs[0, i].set_title(filenames[i])

plt.show() 

輪郭の階層クラスタリング

輪郭を重ね合わせるときフーリエ記述子の差が最小になるようにしました。それでも、輪郭が異なればその差は0にはならず、この値の大きさがつまり形の違い度合いと考えることができます。

そこで、この値を用いて階層クラスタリングをしてみました。今回の3枚の写真に写っているニンジン全てを対象に階層クラスタリングをした結果を以下に示します。

オレンジ色の枝、緑色の枝、赤色の枝にはそれぞれ、短小、細長、中くらいのニンジンがクラスタリングされているように思います。 赤色の枝のクラスタはさらに円筒形または円錐形といった違いでクラスタができているように見えます。

クラスタリングができたからといって嬉しいことは特に思いつきませんが個人的にはおもしろいと思います。

contours = [cnt for _contours in contours_list for cnt in _contours]
indv_imgs = [indv for _indv_imgs in indv_imgs_list for indv in _indv_imgs]

# 輪郭間の距離（どれだけ似てないか、大きい方が似てない）を計算
dist_array = []
n_contours = len(contours)
for i in range(n_contours-1):
    for j in range(i+1, n_contours):
        _, dist = fit.estimateTransformation(contours[j], contours[i])
        dist_array.append(dist)

    
fig = plt.figure()
gs = gridspec.GridSpec(2, n_contours, wspace=0.1, height_ratios=[3, 1])
ax = fig.add_subplot(gs[0, :])
result = linkage(dist_array, method='ward')
d = dendrogram(result, ax=ax)
for i, l in enumerate(d['leaves']):
    ax = fig.add_subplot(gs[1, i])
    ax.imshow(indv_imgs[l])
    ax.axis('off')
    if i == 0:
        ax.tick_params(axis='x', rotation=55)
fig.align_labels()  
plt.show()

フーリエ記述子の主成分分析

フーリエ記述子そのままの値では輪郭の各周波数の成分の強度を表しているとは言えますがあまり解釈しやすい値とは言えないと思います。

そこでよく行われるのが主成分分析です。ここでも主成分分析を行い、第一主成分と第二主成分についてのプロットを書いてみました。

結果は第一主成分で69.5%、第二主成分で15.7%の形のバリエーションを説明できるようです。 また、第一主成分が小さいと細長型、大きいと短小型になるように見えます。 さらに、第二主成分は大きいと円筒型、小さいと円錐型になるように見えます。

fourier_descriptors = np.array([
    cv2.ximgproc.fourierDescriptor(cv2.ximgproc.contourSampling(tcnt, 1024))[1:].ravel()
    for _contours in transformed_contours_list for tcnt in _contours
])
offset_images = [
    OffsetImage(img, zoom=.05)
    for _indv_imgs in indv_imgs_list for img in _indv_imgs
]

pca = PCA(n_components=2)
result = pca.fit_transform(fourier_descriptors)
print(pca.explained_variance_ratio_)

fig, ax = plt.subplots()
artists = []
x, y = result[:, 0], result[:, 1]
for x0, y0, im in zip(x, y, offset_images):
    ab = AnnotationBbox(im, (x0, y0), xycoords='data', frameon=False)
    artists.append(ax.add_artist(ab))
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.update_datalim(np.column_stack([x, y]))
ax.autoscale()
plt.show()

[0.69501555 0.15736549]

輪郭（形）の数値化はなんとなく難しそうなイメージがあったのですが、データ自体は意外と簡単に得られることがわかりました。

参考

• https://pytorch.org/vision/stable/auto_examples/plot_visualization_utils.html#instance-seg-output
• https://docs.opencv.org/4.x/dd/ddc/group__ximgproc__fourier.html