輪郭の曲率を計算してみた

# モジュールのインポート
import cv2
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from matplotlib.animation import FuncAnimation

入力画像

今回はキュウリ画像を例にその輪郭の曲率を計算してみました。

img = cv2.imread('cucumber.jpeg')
fig, ax = plt.subplots()
ax.imshow(img[:, :, ::-1]) # BGR -> RGB
plt.show()

二値化から輪郭抽出まで

まずは二値化し輪郭抽出を行います。一応モルフォロジー演算（Open, Close）をしていますが、結構輪郭がギザギザしていることがわかります。

# Lab色空間のbチャンネルを使ってニ値画像を作成
_, _, b = cv2.split(cv2.cvtColor(img, cv2.COLOR_BGR2Lab))
b = cv2.GaussianBlur(b, (5, 5), 10)
_, mask = cv2.threshold(b, 140, 255, cv2.THRESH_BINARY)

# モルフォロジー演算によるマスクの前処理
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (85, 85))
mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, kernel, iterations=1)
mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel, iterations=1)

# 輪郭抽出
contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
nbElt = 1024
# 等間隔にサンプリング
contour = cv2.ximgproc.contourSampling(contours[0], nbElt=nbElt)

# 輪郭の描画
def plot_contour(img, contour, ax):
    x, y = contour[:, 0, 0], contour[:, 0, 1]
    for a in ax:
        a.imshow(img[:, :, ::-1])
        a.plot(x, y)
    # top
    ax[1].set_xlim(600, 1300)
    ax[1].set_ylim(1200, 0)
    # middle
    ax[2].set_xlim(800, 1500)
    ax[2].set_ylim(2400, 1200)
    # bottom
    ax[3].set_xlim(800, 1500)
    ax[3].set_ylim(3600, 2400)

fig, ax = plt.subplots(1, 4, figsize=(9, 3))
plot_contour(img, contour, ax)
plt.tight_layout()
plt.show()

フーリエ変換による輪郭の平滑化

輪郭を周期関数と見て離散フーリエ変換をしたのち、高周波成分を０にして逆変換することで輪郭を平滑化します。nは使用する係数の数で少ないほど強い平滑化になります。

今回はn=20としました。

def smooth_fourier(contour, n):
    assert n > 0
    dft = cv2.dft(contour)
    dft[n+1:-n] = 0.0
    contour = cv2.idft(dft, flags=cv2.DFT_SCALE)
    contour = cv2.ximgproc.contourSampling(contour, nbElt=len(contour))
    return contour

fig, ax = plt.subplots(3, 4, figsize=(10, 10))
for i, n in enumerate([5, 20, 50]):
    contour_smoothed = smooth_fourier(contour, n)
    plot_contour(img, contour_smoothed, ax[i, :])
    for a in ax[i, :]:
        a.set_title(f'n={n}')
plt.show()

ガウシアンフィルタによる平滑化

一応、他の平滑化手法としてガウシアンフィルタも適用してみました。こちらはsigmaの値が大きいほど強い平滑化になります。

今回はsigma=8としました。

from scipy.ndimage import gaussian_filter1d

def smooth_gaussian(contour, sigma):
    assert sigma > 0
    x = contour[:, 0, 0]
    y = contour[:, 0, 1]
    x = gaussian_filter1d(x, sigma, mode='wrap')
    y = gaussian_filter1d(y, sigma, mode='wrap')
    contour = np.stack([x, y], axis=-1)
    contour = np.expand_dims(contour, 1)

    return contour

fig, ax = plt.subplots(3, 4, figsize=(10, 10))
for i, sigma in enumerate([16, 8, 1]):
    contour_smoothed = smooth_gaussian(contour, sigma)
    plot_contour(img, contour_smoothed, ax[i, :])
    for a in ax[i, :]:
        a.set_title(f'sigma={sigma}')
plt.show()

曲率の計算

いよいよ輪郭の曲率を計算します。輪郭を $\gamma(t) = (x(t), y(t))$ として媒介変数表示したとき、符号付き曲率は次のように定義されます。

$k = \frac{x'y'' - x''y'}{(x'^2 + y'^2)^{3/2}}$

https://en.wikipedia.org/wiki/Curvature より

def calc_curvature(points):
    n = len(points)
    x = points[:, 0, 0]
    y = points[:, 0, 1]

    # 境界部分の微分の安定化のため
    x = np.concatenate([x]*3)
    y = np.concatenate([y]*3)
    
    # 微分
    dx = np.gradient(x)
    dy = np.gradient(y)
    ddx = np.gradient(dx)
    ddy = np.gradient(dy)

    # 曲率
    curvature = (ddx * dy - dx * ddy) / (dx **2 + dy **2)**1.5
    return curvature[n:2*n]

オリジナルの輪郭と、平滑化された輪郭の曲率を計算してみます。オリジナルの輪郭に対する曲率はノイズだらけですが、平滑化することでなんとなくの傾向が見えてきました。

fig, ax = plt.subplots(figsize=(8, 4))

curvature = calc_curvature(contour)
ax.plot(np.roll(curvature, -100), label='original', alpha=0.5)

curvature = calc_curvature(smooth_fourier(contour, n=20))
ax.plot(np.roll(curvature, -100), label='fourier')

curvature = calc_curvature(smooth_gaussian(contour, sigma=8))
ax.plot(np.roll(curvature, -100), label='gaussian')

ax.legend()
plt.show()

続いて輪郭の各位置と曲率の対応関係を可視化するため、動画として出力してみました。輪郭と曲率のグラフについて、それぞれ同じインデックスの点を動く点として表示しています。

def plot_curvature(img, contour, curvature):
    fig, (ax_ct, ax_cv) = plt.subplots(1, 2)

    # 入力画像と輪郭
    x, y = contour[:, 0, :].T
    ax_ct.imshow(img[:, :, ::-1])
    ax_ct.plot(x, y, 'tab:blue')

    n = len(contour)
    m = n // 2
    i = np.arange(n)
    # 輪郭上の動点
    point_ct, = ax_ct.plot(x[m], y[m], 'tab:blue', marker='o')
    # 曲率
    line_cv, = ax_cv.plot(i, curvature, 'tab:blue')
    # 曲率上の動点
    point_cv, = ax_cv.plot([m], [curvature[m]], 'tab:blue', marker='o')

    def update(frame):
        frame = frame
        rolled = np.roll(curvature, -frame)
        point_ct.set_data([x[(m+frame)%n]], [y[(m+frame)%n]])
        line_cv.set_data(i, rolled)
        point_cv.set_data([m], [rolled[m]])
        return point_ct, line_cv, point_cv,

    ani = FuncAnimation(fig, update, frames=range(0, n, 10), interval=5, blit=True)
    return ani

contour = smooth_fourier(contour, 20)
curvature = calc_curvature(contour)
ani = plot_curvature(img, contour, curvature)
ani.save('curvature.gif', writer='imagemagick', fps=10)

こう見ると曲率が極大を取るところはキュウリのヘタと先端部分であることがわかります。このことから、キュウリのヘタと先端部分を抽出し、プロットすることも可能となります。

def where_tip(contour):
    curvature = calc_curvature(contour)

    # 極大値が輪郭の配列先頭にあると極大として見えないのでずらす
    min_i = curvature.argmin()
    curvature_rolled = np.roll(curvature, -min_i)
    contour_rolled = np.roll(contour, -min_i, axis=0)

    peaks, _ = find_peaks(curvature_rolled)
    assert len(peaks) >= 2
    peaks = sorted(peaks, key=lambda p: curvature_rolled[p], reverse=True)[:2]

    # 曲率極大を示す輪郭位置の座標
    tips = contour_rolled[peaks]

    return tips

tip1, tip2 = where_tip(contour)

fig, ax = plt.subplots()
ax.imshow(img[:, :, ::-1])
ax.plot(contour[:, :, 0], contour[:, 0, 1])
ax.plot(tip1[0, 0], tip1[0, 1], 'r*')
ax.plot(tip2[0, 0], tip2[0, 1], 'c*')
plt.show()