OpenCVを利用した葉面積の推定

• このページのJupyter Notebookファイル
• google colaboratoryで開く
• 画像データ
• 解説動画

Pythonで画像を扱う際に、OpenCVというライブラリがよく利用されます。 今回はこのOpenCVを用いて、作物成長段階の指標となる葉面積を推定します。

まずはOpenCVモジュールを読み込みます。 以下のチャンクを実行すると、今動いているPythonの環境にOpenCVの機能が追加されます。

import cv2

作業しているパソコンにOpenCVがインストールされていない場合は、エラーが出ますので、以下のコマンドでインストールして、再度実行してください。

!pip install opencv-python

Requirement already satisfied: opencv-python in /home/takahisa/anaconda3/envs/py38/lib/python3.8/site-packages (4.4.0.44)
Requirement already satisfied: numpy>=1.17.3 in /home/takahisa/anaconda3/envs/py38/lib/python3.8/site-packages (from opencv-python) (1.19.1)

うまくインストールされれば、以下のコマンドでバージョンが表示されます。

cv2.__version__

'4.4.0'

OpenCVでの画像の表示

まずは画像ファイルをJupyter Notebookに読み込んで表示してみましょう。

画像は、.ipynbファイルと同じフォルダにdataという名前のフォルダを作成し、その中に保存してください。以下では、www.pxfuel.comから取得した画像を利用しています。

# ファイルへのパスをimage_path変数に代入する。
image_path = 'data/field.jpg'

# 画層をOpenCVのimage形式で読み込む
img = cv2.imread(image_path)

次のチャンクを実行すると、imageという名前の別ウィンドウが立ち上がり、画像が表示されます。(Jupyter labをWSL2上やGoogle Colaboratoryで動かしている場合、下のコードは機能しないかもしれません。参考:Can’t use X-Server in WSL 2)

表示されたウィンドウは、キーボードのいづれかのキーを押して閉じてください。

# OpenCV のimshow関数で、画像を表示する
# 作業環境によっては動かない
# 動かない場合は飛ばしてよい
cv2.imshow('image',img)
cv2.waitKey(0)
cv2.destroyAllWindows()

うまく実行されれば別ウィンドウで次の画像が表示されます。

表示されたウィンドウを閉じてもセルの番号が[*]と表示されたままの場合は、Jupyter NotebookメニューのKernel -> Restart Kernelで、カーネルを再起動してください。

画像の表示はこの方法でもいいですが、matplotlibというライブラリを利用して、Jupyter Notebook内にインラインで表示させる下記の方法がよく用いられます。

まずはmatplotlibをインポートします。(必要であればmatplotlibを!pip install matplotlib からインストールしてください。)

import matplotlib.pyplot as plt

ここでは、matplotlibライブラリ内のpyployというサブモジュールをpltという別名でインポートしました。このように別名を使ってインポートすることにより、たとえばmatplotlib.pyplot.imshow()と書く代わりに、plt.imshow()と書くことができるようになります。

import文はプログラムの一番最初にまとめて記述することが慣例ですが、ここでは必要になる箇所で適宜宣言することにします。

では、matplotlibのmatshow関数を使って画像を表示しましょう。

plt.matshow(img)

<matplotlib.image.AxesImage at 0x7f66c35541f0>

色が変ですね。

これは、OpenCVとmatplotlibでのデフォルトのRGBチャンネルの順序が異なるため、赤と青が入れ替わってしまっているからです。

OpenCVはBGRの順、matplotlibはRGBの順でカラーを取り扱います。 正しく読み込むには以下のように読み込みます。

# 画層をRGBの順で読み込む
img_bgr = cv2.imread(image_path)

# cvtColor() 関数でBGRからRGBへ変換する
img_rgb = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2RGB)
# matplotlibで表示
plt.matshow(img_rgb)

<matplotlib.image.AxesImage at 0x7f66b8443790>

.jpgファイルとして保存された写真をPythonで表示することができました。

画像データとはなにか

写真や衛星画像のようなデータはラスタデータと呼ばれ、一定の間隔で格子状に区切られた各グリッドに値をもったデータ形式です。それぞれのグリッドのことをピクセルと呼ぶので、ピクセルデータと呼ぶこともあります。上で画像をプロットした際に、X軸とY軸に数値が振られていますが、この数値が画像のピクセルを表しています。通常左上を起点とした座標で表されます。

ラスタデータと対になるデータはベクタデータです。分野によっては、ベクタデータは、ベクトル、ポリゴンなどとも呼ばれます。

詳しくはESRIのサイトに解説があるので、参考にしてください。

• ベクターデータとは: www.esrij.com

• ラスターデータとは: wwww.esrij.com

今回読み込んだ画像はラスタデータです。ここではPythonの中でどのように保存されているかを見ていきます。

まずはcv2.imread()で読み込んだimg_bgrをprint()関数を用いて表示させてみます。

print(img_bgr)

[[[ 71 205 152]
  [ 68 202 149]
  [ 68 199 148]
  ...
  [ 43  66  92]
  [ 68  95 121]
  [105 136 161]]

 [[ 70 206 154]
  [ 68 204 152]
  [ 64 200 148]
  ...
  [ 58  85 111]
  [ 78 109 134]
  [107 143 167]]

 [[ 67 206 155]
  [ 64 203 152]
  [ 59 198 147]
  ...
  [ 76 113 135]
  [ 85 126 149]
  [102 146 169]]

 ...

 [[187 217 252]
  [191 219 254]
  [181 206 246]
  ...
  [104 120 143]
  [ 98 103 128]
  [ 90  91 117]]

 [[168 197 224]
  [184 210 240]
  [189 214 246]
  ...
  [115 136 157]
  [104 115 137]
  [ 90  96 119]]

 [[150 175 195]
  [163 187 209]
  [167 190 216]
  ...
  [117 142 158]
  [101 116 135]
  [ 82  93 113]]]

print()関数によって、画像の中のデータが表示されました。画像データは、整数の値が入ったリストが入れ子になった構造をしていることがわかります。

type()関数で、データ型を確認します。

type(img_bgr)

numpy.ndarray

img_bgrの実態は、numpy.ndarrayクラスのデータであることが分かりました。 numpy.ndarrayクラスは、shapeプロパティを持つので、これを確認してみます。

img_bgr.shape

(853, 1280, 3)

実は、OpenCVで読みこんんだデータは、各ピクセルにおける輝度が[青、緑、赤]の順で保存され、その各ピクセルの値は、画像の左上から右上、さらに上から下の順で3重の入れ子構造で保存されています。そのため、shapeを確認することで、読み込んだ画像は高さ853ピクセル、幅1280ピクセルのサイズで、3チャンネルの画像であったことがわかりました。

今回の画像の例では、print()関数で一行目に表示された[ 71 205 152]の値は、左上のピクセルの青色の値が71、緑色が205、赤色が152の値を持っています。値が大きいほど、ディスプレイ上でその色が明るく光るので、このピクセルは緑>赤>青の順の強さで光って、黄緑色であることが推測できます。

cv2.cvtColor()関数で変換されたimg_rgbのデータと比べてみましょう。

print(img_rgb)

[[[152 205  71]
  [149 202  68]
  [148 199  68]
  ...
  [ 92  66  43]
  [121  95  68]
  [161 136 105]]

 [[154 206  70]
  [152 204  68]
  [148 200  64]
  ...
  [111  85  58]
  [134 109  78]
  [167 143 107]]

 [[155 206  67]
  [152 203  64]
  [147 198  59]
  ...
  [135 113  76]
  [149 126  85]
  [169 146 102]]

 ...

 [[252 217 187]
  [254 219 191]
  [246 206 181]
  ...
  [143 120 104]
  [128 103  98]
  [117  91  90]]

 [[224 197 168]
  [240 210 184]
  [246 214 189]
  ...
  [157 136 115]
  [137 115 104]
  [119  96  90]]

 [[195 175 150]
  [209 187 163]
  [216 190 167]
  ...
  [158 142 117]
  [135 116 101]
  [113  93  82]]]

cv2.cvtColor()によって、赤色の値と青色の値が入れ替えられたことが確認できます。

また、numpy.arrayは以下のようにすることで、データの一部を取り出すことができます。

array[1軸目の範囲, 2軸目の範囲, 3軸目の範囲, ...]

今回の画像データは縦、横、色チャンネルの3次元データです。左上の縦100x横200ピクセルに対応するimg_rgbは次のように取り出すことができます。

img_rgb[0:100, 0:200, :]  # 3軸目の範囲":" は"すべての範囲"を意味する

array([[[152, 205,  71],
        [149, 202,  68],
        [148, 199,  68],
        ...,
        [125, 144,  89],
        [127, 150,  96],
        [128, 154, 106]],

       [[154, 206,  70],
        [152, 204,  68],
        [148, 200,  64],
        ...,
        [142, 160, 108],
        [143, 163, 112],
        [139, 162, 116]],

       [[155, 206,  67],
        [152, 203,  64],
        [147, 198,  59],
        ...,
        [152, 168, 123],
        [146, 163, 119],
        [135, 156, 113]],

       ...,

       [[147, 187, 101],
        [137, 176,  97],
        [137, 176, 109],
        ...,
        [145, 194,  42],
        [154, 204,  47],
        [168, 219,  54]],

       [[147, 188,  96],
        [142, 182,  96],
        [143, 181, 108],
        ...,
        [150, 199,  45],
        [155, 205,  46],
        [163, 216,  50]],

       [[149, 190,  98],
        [145, 185,  99],
        [143, 180, 110],
        ...,
        [150, 199,  48],
        [150, 201,  46],
        [156, 208,  48]]], dtype=uint8)

この部分だけをプロットしてみましょう。

plt.imshow(img_rgb[0:100, 0:200, :])

<matplotlib.image.AxesImage at 0x7f66b83b4610>

画像の左上あたりが切り取られていることが確認できました。

チャンネルの分解

OpenCVにあるsplit()関数を用いると、上記の画像形式のnumpy.ndarrayをチャンネルごとに分割することができます。

r, g, b = cv2.split(img_rgb)

cv2.spolit()関数は結果を長さ3のタプルで返します。返り値がタプルの場合、代入記号=の左側にタプルの要素と同じ数の変数を置くことによって、タプルの各要素をそれぞれの変数に代入することができます。これをアンパックといいます。参考:Pythonの関数で複数の戻り値を返す方法

cv2.spolit(img_rgb)は、一つ目の返り値に赤チャンネル、二つ目に緑チャンネル、三つ目に青チャンネルを返すので、それぞれのチャンネルの値をr、g、bに代入しました。

各チャンネルのtype、shapeを確認しましょう。

type(r)

numpy.ndarray

r.shape

(853, 1280)

print(r)

[[152 149 148 ...  92 121 161]
 [154 152 148 ... 111 134 167]
 [155 152 147 ... 135 149 169]
 ...
 [252 254 246 ... 143 128 117]
 [224 240 246 ... 157 137 119]
 [195 209 216 ... 158 135 113]]

元のimg_rgbは3次元のnumpy.ndarrayでしたが、各チャンネルは画像の縦と横の2次元のnumpy.ndarray形式になったことが確認できました。

matplotlib.pyplotを用いて、チャンネルごとに画像としてプロットしてみましょう。

plt.matshow(r, cmap='Reds')  # 赤色の画像で表示する
plt.colorbar()               # 画像の右側にスケールバーを表示する
plt.show()

元の画像と各チャンネルの画像を並べてプロットするには以下のように行います。

channels = [img_rgb, r, g, b]
colors = (None, 'Reds','Greens','Blues')
titles = ('Original', 'Red channel', 'Green channel', 'Blue channel')

fig, axs = plt.subplots(1, 4, figsize=(12,3))  # figとaxを作成

# axs[0].set_title('Original')
# axs[1].set_title('Red channel')
# axs[2].set_title('Grenn channel')
# axs[3].set_title('Blue channel')

# axs[0].matshow(img_rgb)
# axs[1].matshow(r, cmap='Reds')
# axs[2].matshow(g, cmap='Greens')
# axs[3].matshow(b, cmap='Blues')

# axs[0].xaxis.tick_bottom()
# axs[1].xaxis.tick_bottom()
# axs[2].xaxis.tick_bottom()
# axs[3].xaxis.tick_bottom()

# forループで各画像を描写する
for idx, (title, channel, color) in enumerate(zip(titles, channels, colors)):
    axs[idx].set_title(title)
    ax = axs[idx].matshow(channel, cmap=color)
    axs[idx].xaxis.tick_bottom()
    if color:
        fig.colorbar(ax, ax=axs[idx], fraction=0.04)

# プロット(axs)間のスペースを調整する
plt.tight_layout()

plt.show()

HSVへの変換

色を表現する方法は、赤・緑・青の各バンドの強さを示すほかに、色相(Hue)・彩度(Saturation)・明度(Value)の3つの値で示す方法があります。これをRGB形式に対して、HSV形式といいます。

画像解析を行う際に、RGBからHSVに変換することがよく行われます。RGBはディスプレイ上で色を表現するには都合がよいですが、各ピクセルが「赤っぽい」のかや、「鮮やかな色」なのかといったことがわかりにくい表現形式になっています。今回のように、画像から「緑色」のピクセルを探す場合、RGBの形式で単純にGの値が大きいところを選択すると、白っぽいところ、つまりRもGもBも大きな値のピクセルも選択されてしまいます。逆に緑っぽいピクセルでも陰になっているとGの値が小さくなるので選択することができません。

HSVに変換することで、色相による選択を容易にすることができるようになります。

RGB形式のデータをHSV形式のデータに変換するには、cv2.cvtColor()関数を利用します。

img_hsv = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2HSV)
img_hsv.shape

(853, 1280, 3)

img_hsvも同様にhとsとvに分割することができます。

h, s, v = cv2.split(img_hsv)

作成したHSVデータをプロットしてみましょう。

channels = [img_rgb, h, s, v]
colors = (None, 'hsv','Greys','Greys')
titles = ('Original', 'Hue', 'Saturation', 'Value')

fig, axs = plt.subplots(1, 4, figsize=(12,3))  # figとaxを作成

# axs[0].set_title('Original')
# axs[1].set_title('Hue')
# axs[2].set_title('Saturation')
# axs[3].set_title('Value')

# axs[0].matshow(img_rgb)
# axs[1].matshow(h, cmap='hsv')
# axs[2].matshow(s, cmap='Greys')
# axs[3].matshow(v, cmap='Greys')

# axs[0].xaxis.tick_bottom()
# axs[1].xaxis.tick_bottom()
# axs[2].xaxis.tick_bottom()
# axs[3].xaxis.tick_bottom()

# forループで各画像を描写する
for idx, (title, channel, color) in enumerate(zip(titles, channels, colors)):
    axs[idx].set_title(title)
    ax = axs[idx].matshow(channel, cmap=color)
    axs[idx].xaxis.tick_bottom()
    if color:
        fig.colorbar(ax, ax=axs[idx], fraction=0.04)

plt.tight_layout()
plt.show()

閾値による分類

画像をHSVへ変換することができたので、色相(Hue)の値が特定の範囲にあるピクセルを抽出することで、葉っぱが写っている面積を算出します。

今回はhの値が40から60の値を抽出してみましょう。

h_min = 40
h_max = 60
img_bin = (h_min < h) & (h < h_max) * 1
plt.matshow(img_bin, cmap='binary_r')
plt.colorbar(ticks=[0, 1])
plt.show()

元の画像で緑色のピクセルを抽出することができています。作成したimg_binが1の値の個所が抽出されたピクセルで、上の画像では白く表示されています。

抽出されたピクセルを、画像のすべてのピクセル数で割ることで、画像に占めるダイズの葉の割合を計算することができます。

print(f'Mask ratio is {img_bin.sum()/img_bin.size}')

Mask ratio is 0.2845270369284877

実測した葉面積とこの画像に対するダイズの葉の割合との校正曲線を作成することで、定点観測により撮影した画像から葉面積を推定することができます。

おまけ

HSVへの画像変換を行わない場合、画像の抽出がどのようになるか確認してみましょう。

g_min = 150
g_max = 250
img_rgb_bin = (g_min < g) & (g < g_max) *1
# plt.imshow(img_bin, cmap='binary_r')
plt.matshow(img_rgb_bin, cmap='binary_r')
plt.colorbar(ticks=[0, 1])
plt.show()

g_minの値を様々に変えてみても、HSV画像で抽出した方が少ないノイズで抽出できています。

また以下のグラフで、HSVで抽出されたピクセルはRGBの値でどのように表現されているものかを確認することができます。

Gの値のみで抽出することが難しいことがわかりますね。

import numpy as np
import random
def plot_3d(img, sampling=1000, mask=None, **kwargs):
    xlab, ylab, zlab = 'R', 'G', 'B'
    y,x,z = img.shape
    img_flat = np.reshape(img, (x*y,z))
    if mask is not None:
        mask_flat = np.reshape(mask, (x*y)).astype('bool')
        img_flat = img_flat[mask_flat, :]

    if sampling != 0:
        # random sample
        img_flat = img_flat[np.random.choice(len(img_flat), sampling, replace=False)]
    colors = img_flat / 255
    fig = plt.figure(**kwargs)
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(xs=img_flat[:,0],
               ys=img_flat[:,1],
               zs=img_flat[:,2],
               s=10,
               c=colors,
               lw=0)
    ax.set_xlabel(xlab)
    ax.set_ylabel(ylab)
    ax.set_zlabel(zlab)
    plt.show()

# matplotlibのインタラクティブプロットを有効化
# うまく行かない場合はコメントアウトしてください。
%matplotlib widget

plot_3d(img = img_rgb,
        mask=img_bin,
        sampling=1000,
        figsize=(10,10))

# インタラクティブを無効化
%matplotlib inline