[DepthAnything][深度推定]ドローン画像の樹頂点と樹冠の検出

はじめに

ドローンや航空機から撮影ができるオルソ画像をもとに木の樹頂点と樹冠の検出をご紹介します。

一般的な森林の解析では、オルソ画像と合わせてレーザー計測によって作成したDSM（撮影物の表面の高さ情報）とDEM（地形の標高モデル）の差分からCHM（撮影物の高さモデル）を作成した上で、樹頂点や樹冠の検出を実行します。

今回は、Depth Anythingを使用することで、オルソ画像のみで樹頂点検出と樹冠分離をやってみたいと思います。

目次

1. DepthAnythingとは
2. 入力データ
3. 深度推定
4. 樹頂点の検出
5. 樹冠の分離
6. おわりに

DepthAnythingとは

Depth Anythingとは深度推定モデルであり、画像から奥行きを推定することのできるモデルです。
2024年6月にV2がリリースされています。
https://github.com/DepthAnything/Depth-Anything-V2

今回はこちらを利用して、オルソ画像からCHMを作成します。

入力データ

DeepForestさんよりサンプルデータ（オルソ画像）をダウンロードし検証いたします。
https://www.df-webservice.com/download

そのままでは、深度推定できないので、rasterioを使用してTiffをPNG変換します。

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from skimage import io, filters, feature
import geopandas as gpd
from shapely.geometry import Polygon, Point
from skimage.measure import label, regionprops, find_contours, approximate_polygon
from skimage.segmentation import watershed

# ファイルパス
input_file = '/data/Ortho.tif'
output_file = 'data/Ortho_RGBA.png'
# 16ビット画像データを8ビットにスケーリングする関数
def scale_to_8bit(img, max_val=None):
    if max_val is None:
        max_val = img.max()
    img = (img / max_val) * 255
    return img.astype(np.uint8)

# GeoTIFFファイルを開いて処理
with rasterio.open(input_file) as src:
    # バンド1, 2, 3をそれぞれ赤、緑、青として読み込む
    red = src.read(1)
    green = src.read(2)
    blue = src.read(3)
    
    # 各バンドを8ビットにスケーリング
    red_8bit = scale_to_8bit(red)
    green_8bit = scale_to_8bit(green)
    blue_8bit = scale_to_8bit(blue)
    
    # アルファチャンネルを作成（不透過の255で設定）
    alpha = np.full(red.shape, 255, dtype=np.uint8)

    # RGBA画像を作成
    rgba = np.dstack((red_8bit, green_8bit, blue_8bit, alpha))

    # RGBA画像をPNGとして保存
    plt.imsave(output_file, rgba)

深度推定

Depth Anythingを使用して深度推定します。

import cv2
import torch

from depth_anything_v2.dpt import DepthAnythingV2

DEVICE = 'cuda' if torch.cuda.is_available() else 'mps' if torch.backends.mps.is_available() else 'cpu'

model_configs = {
    'vits': {'encoder': 'vits', 'features': 64, 'out_channels': [48, 96, 192, 384]},
    'vitb': {'encoder': 'vitb', 'features': 128, 'out_channels': [96, 192, 384, 768]},
    'vitl': {'encoder': 'vitl', 'features': 256, 'out_channels': [256, 512, 1024, 1024]},
    'vitg': {'encoder': 'vitg', 'features': 384, 'out_channels': [1536, 1536, 1536, 1536]}
}

encoder = 'vitl' # or 'vits', 'vitb', 'vitg'

model = DepthAnythingV2(**model_configs[encoder])
model.load_state_dict(torch.load(f'checkpoints/depth_anything_v2_{encoder}.pth', map_location='cpu'))
model = model.to(DEVICE).eval()

raw_img = cv2.imread('your/image/path')
depth = model.infer_image(raw_img) # HxW raw depth map in numpy

元のオルソ画像

深度推定結果

なかなかきれいに推定できたのではないでしょうか。

樹頂点の検出

画像の読み込みと前処理を実行します。

# 0. RAW画像の読み込み
image_path = './Ortho_RGBA.png'
raw_image = io.imread(image_path)

# 1. CHM画像の読み込み
image_path = './chm.png'
chm_image = io.imread(image_path)

# 2. 高さによるフィルタリング
height_min = 100
height_max = 255
filtered_chm = np.where((chm_image >= height_min) & (chm_image <= height_max), chm_image, 0)

# 2. 画像の平滑化 (ノイズを減らすためにガウスフィルタを適用)
smoothed_chm = filters.gaussian(filtered_chm, sigma=3)

樹頂点を検出します。

# 3. 樹頂点の検出 (局所的な極大値を検出)
local_maxi = feature.peak_local_max(smoothed_chm, min_distance=20, threshold_rel=0.1)

points = [Point(xy[1], xy[0]) for xy in local_maxi]  # (行, 列)を(x, y)に変換
gdf_treetop = gpd.GeoDataFrame(geometry=points)

fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(chm_image, cmap='gray')
ax.plot(local_maxi[:, 1], local_maxi[:, 0], 'r.', markersize=5)  # Tree tops as red dots

plt.show()

min_distance（検出した点同士の最低距離）とthreshold_relは、調整が必要になります。

樹冠の分離

watershed segmentationを使用して、樹頂点を中心に領域をセグメンテーションします。

# 4. Watershedアルゴリズムを使ったセグメンテーション
# 樹頂点をマーカーとしてWatershedを適用
markers = np.zeros_like(smoothed_chm, dtype=np.int32)
markers[local_maxi[:, 0], local_maxi[:, 1]] = np.arange(1, len(local_maxi) + 1, dtype=np.int32)


# 画像の勾配（セグメンテーション用のランドスケープ）を計算
gradient = filters.sobel(smoothed_chm)

# Watershedセグメンテーションの適用
labels = watershed(gradient, markers, mask=(filtered_chm > 0))

# Plot
fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(chm_image, cmap='gray')
ax.imshow(labels, cmap='nipy_spectral', alpha=1)  # Display the segmented regions with transparency
ax.plot(local_maxi[:, 1], local_maxi[:, 0], 'r.', markersize=5)  # Tree tops as red dots

plt.show()

仕上げに抽出した領域をシェープファイルに変換します。

# 5. ポリゴンデータの生成
# ラベル付けされた画像から各領域を抽出
label_image = label(labels)
padded_labels = np.pad(label_image, pad_width=1, mode='constant', constant_values=0)

# 有効なポリゴンを保存するリスト
valid_polygons = []

# 各領域の輪郭を抽出し、ポリゴンを作成
for region in regionprops(label_image):
    if region.area > 10:  # ある程度大きな領域のみを対象とする
        contours = find_contours(padded_labels == region.label, 0.5)
        for contour in contours:
            # 輪郭を単純化してポリゴンに変換
            simplified_contour = approximate_polygon(contour, tolerance=2.5)
            if len(simplified_contour) >= 3:  # ポリゴンとして最低3点が必要
                poly = Polygon(simplified_contour[:, [1, 0]])  # xとyを入れ替えてポリゴンを作成
                valid_polygons.append(poly)

# 6. GeoDataFrameの作成
# 有効なポリゴンからGeoDataFrameを作成
gdf_crown = gpd.GeoDataFrame(geometry=valid_polygons)

# Plot the resulting shapefile
fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(raw_image)
gdf_crown.boundary.plot(ax=ax, color='blue')
gdf_treetop.plot(ax=ax, color='red')
plt.show()

それっぽい結果を出力することができました。

おわりに

オルソ画像からだけでも、精度をそこまで気にしなければ、樹頂点検出と樹冠分離ができることを確認できました。

もっと精緻な解析が必要な場合は、サンプルデータの提供元であるDeepForestさんのDF Scannerをご利用されるのが良いかと思います。

参考