輪郭抽出

import urllib.request

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import rgb_to_hsv
from scipy import ndimage as ndi

base = "https://assets.ku-compbio.noshita.net/2026/P2-04/images/"
filename = "Pedu1_1_4.png"
urllib.request.urlretrieve(base + filename, filename)

img = plt.imread(filename)[:, :, :3]
img_gray = img @ [0.2989, 0.587, 0.114]
plt.imshow(img)

plt.hist(img_gray.ravel(), bins=256)

img_bin = img_gray < 0.7
plt.imshow(img_bin, cmap="gray")

def otsu_threshold(gray, nbins=256):
    """輝度ヒストグラムからクラス間分散が最大となる閾値を返す."""
    hist, edges = np.histogram(gray.ravel(), bins=nbins, range=(0.0, 1.0))
    centers = (edges[:-1] + edges[1:]) / 2
    weight0 = np.cumsum(hist)  # 閾値以下の画素数
    weight1 = weight0[-1] - weight0  # 閾値より上の画素数
    sum0 = np.cumsum(hist * centers)
    mean0 = np.divide(sum0, weight0, out=np.zeros(nbins), where=weight0 > 0)
    mean1 = np.divide(sum0[-1] - sum0, weight1, out=np.zeros(nbins), where=weight1 > 0)
    var_between = weight0 * weight1 * (mean0 - mean1) ** 2
    return centers[np.argmax(var_between)]


t = otsu_threshold(img_gray)
t

plt.imshow(img_gray < t, cmap="gray")

hsv = rgb_to_hsv(img)
hue = hsv[:, :, 0]
sat = hsv[:, :, 1]

# 緑色相（およそ 0.05〜0.55）かつ彩度のある画素．厳しくすると葉の基部が削れるため広めにとる
green = (hue > 0.05) & (hue < 0.55) & (sat > 0.06)
plt.imshow(green, cmap="gray")

leaf = ndi.binary_opening(green)  # 細かなノイズを除去（強すぎると葉の基部が削れる）
leaf = ndi.binary_fill_holes(leaf)  # 内部の穴を埋める
plt.imshow(leaf, cmap="gray")

fig, ax = plt.subplots()
ax.imshow(img)
cs = ax.contour(leaf, levels=[0.5], colors="red", linewidths=1.5)

contours = cs.allsegs[0]  # 各輪郭は (N, 2) の (x, y) 座標
len(contours)

def shoelace_area(xy):
    """順序付き座標 (N, 2) が囲む多角形の面積（Shoelace formula）."""
    x, y = xy[:, 0], xy[:, 1]
    return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))


min_area = 0.005 * leaf.size
leaf_contours = [c for c in contours if shoelace_area(c) > min_area]
[int(shoelace_area(c)) for c in leaf_contours]

fig, ax = plt.subplots()
ax.imshow(img)
for c in leaf_contours:
    ax.plot(c[:, 0], c[:, 1], "r-", lw=1.5)

biggest = max(leaf_contours, key=shoelace_area)
plt.plot(biggest[:, 0], -biggest[:, 1])  # 画像座標に合わせて y を反転
plt.axis("equal")