【scikit-image】画像の特定方向のプロファイルを簡単に取得する方法

はじめに

画像処理において、任意の直線上の強度プロファイル（プロジェクション）の作成は、画像の特定の方向に関する特性を理解するためによく使われる方法の一つです。しかし、一般的なプロファイル処理には、領域の切り取りや回転行列の準備、プロファイルの垂直方向の幅など、さまざまな要素が組み合わさるため、自作するにはプログラムがやや複雑であり、バグや計算コストの増加が懸念されます。

そこで、公開ライブラリを利用することを考えますが、scikit-imageにはその処理に特化したprofile_line()が用意されているため、本記事ではこの関数の使い方や補足事項を解説します。

方法

使用したコードはGoogle Colabこちらに掲載していますので、よろしければご覧ください。

モジュールのインポート

本記事で使用するモジュールは以下の通りです。skimage.measureのprofile_lineが本記事で使用するプロファイル作成用の関数です。

使用するモジュールのインポート

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from skimage.measure import profile_line
from scipy.stats import multivariate_normal

本記事で使用したscikit-imageのバージョンは0.19.3です。

サンプル画像の用意

本研究で使用するサンプル画像を準備します。ここでは、多変量ガウス分布にポアソンノイズを加えた画像を使いますが、関数には特に記事の本題とは関連しないので説明は割愛します。

サンプル画像作成用の関数

def get_sample_image():
    y, x = np.mgrid[-50:51:, -100:101:]
    pos = np.dstack((y, x))
    mean = np.array([0, 0])
    sigma_y, sigma_x = 10, 30
    sigma_yx = sigma_y*sigma_x*0.9
    cov = np.array([[sigma_y**2, sigma_yx], [sigma_yx, sigma_x**2]])
    rv = multivariate_normal(mean, cov)
    img = rv.pdf(pos)
    
    # ポアソン分布に従うランダムサンプリング
    np.random.seed(2024)
    num_samples = 3000
    img_samples = np.random.poisson(img * num_samples)
    
    return img_samples

強度プロファイルの作成と描画

プロファイルの作成には、profile_line()を使います。また、描画の際はプロファイル領域と画像の対応がうまく取れるように実装し、matplotlibで描画します。以下がその実装コードになります。

実装コード

# プロファイルの始点と終点と幅
start = np.array([30, 100])  # プロファイルの始点 (y, x)の順
end = np.array([80, 85])  # プロファイルの終点 (y, x)の順
lw = 50  # プロファイルの幅

# 画像を取得
img = get_sample_image()

# プロファイルを取得
profile = profile_line(img, start, end, linewidth=lw, reduce_func=np.mean)
profile_region = profile_line(img, start, end, linewidth=lw, reduce_func=None)

# プロファイル情報を描画
# ベクトルに直交する単位ベクトルを計算
vec = end - start
ortho_vec = np.array([vec[1], -vec[0]]) 
ortho_unit_vec = ortho_vec / np.linalg.norm(ortho_vec)

# プロファイルの幅に関する長方形の頂点の計算
rectangle_points = [start + ortho_unit_vec * lw / 2,
                    start - ortho_unit_vec * lw / 2,
                    end - ortho_unit_vec * lw / 2,
                    end + ortho_unit_vec * lw / 2,
                    start + ortho_unit_vec * lw / 2]

# 座標系を揃えるためyとxの対応を反転
rectangle_points = np.array(rectangle_points)[:, [1, 0]]

# グラフ描画
fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(15, 5))

# 元画像
ax1.imshow(img, origin='lower')
rectangle = patches.Polygon(rectangle_points, closed=True, fill=True, color='r', alpha=0.1)
ax1.add_patch(rectangle)
ax1.quiver(start[1], start[0], vec[1], vec[0], angles='xy', scale_units='xy', scale=1, color='r')
ax1.set_title('original image')
ax1.set_xlabel('x')
ax1.set_ylabel('y')

# プロファイル領域
ax2.set_title('profile region')
ax2.imshow(np.swapaxes(profile_region, 0, 1), origin='lower')
ax2.set_xlabel('distance $\\to$')

# プロファイル結果
ax3.plot(profile)
ax3.set_title('profile result')
ax3.set_xlabel('distance $\\to$')
ax3.set_ylabel('value')

# レイアウト調整
fig.tight_layout()
plt.show()

図の説明

• 左図：元画像とプロファイル領域。矢印がプロファイルの方向、赤枠がプロファイルの幅を表す。
• 中央図：プロファイル領域を表し、横軸が左図の矢印に対応する。
• 右図：プロファイルの結果。横軸が左図の矢印に対応する。

コードの説明

• start（始点）、end（終点）でプロファイルの範囲を(y,x)の順で指定します。endで指定した端の値も含みます。例えば、end = np.array([80, 85])としたとき、y座標の80、x座標の85をデータ点として含みます。この操作はスライスの指定方法(その値を含まない)と異なるので注意が必要です。上図左がそのプロファイルの位置を示しており、(y,x)の順序で、始点(30,100)、終点(80,85)となっています。
• lwでプロファイルの幅を指定します。実装例では、lw=50とし、上図左で赤の四角がその幅に対応します。
• profile = profile_line(img, start, end, linewidth=lw, reduce_func=np.mean)で、プロファイルを作成します。reduce_funcは、lwで指定した幅の計算種類の関数を表し、例えば、reduce_func=np.meanとすると、平均を返します(上図右がその結果)。また、reduce_func=Noneとすると、プロファイルの領域の計算をする前の画像そのものを返します(上図中央がその結果)。

補足

ここでは、profile_line()のいくつか主要な引数について補足します。

reduce_funcについて

reduce_funcには、np.meanの他にもnp.sumやnp.sqrtなど自由に設定できます。例えば、下図の左にnp.mean(平均)、右にnp.sum(合計)としたときの結果を表示します。

profile_mean = profile_line(img, start, end, linewidth=lw, reduce_func=np.mean)
profile_sum = profile_line(img, start, end, linewidth=lw, reduce_func=np.sum)

# グラフ描画
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10, 5))

# プロファイル結果(平均)
ax1.plot(profile_mean)
ax1.set_title('profile mean result')
ax1.set_xlabel('distance $\\to$')
ax1.set_ylabel('value')

# プロファイル結果(合計)
ax2.plot(profile_sum)
ax2.set_title('profile sum result')
ax2.set_xlabel('distance $\\to$')
ax2.set_ylabel('value')

# レイアウト調整
fig.tight_layout()
plt.show()

プロファイル領域が画像をはみ出るとき

profile_line()関数は、プロファイル領域が画像をはみ出る場合に補完処理が行われた上で計算される機能を持ち、この設定はprofile_line()のmode引数から設定できます。

• mode='reflect'（デフォルト）
  
  中央図で、画像のない部分が'refrect'(反射)補完されていることがわかります。右図では、そのような補完が行われたデータのプロファイルになっています。デフォルトがこのモードになっている点に留意してください。

• mode='constant'

mode='constant'とすると、画像のはみ出た部分を特定の値で埋めます。その特定の値の設定は、profile_line()のcval引数で指定します。この例ではcval=0としています。

スプライン補完

profile_line()関数のorder引数から、スプライン補完の次数を設定できます。スプライン補完の説明については、pythonでxy座標上の離散点をスプライン補間などが参考になります。

デフォルトは、order=0(入力画像の型がbool)、order=1(入力画像の型がbool以外)となっています。
orderに指定できる次数は、0, 1, 2, 3, 4, 5となります。以下に、スプライン補完のorderを変えた際の結果を示します。スプライン補完の次数が大きほど得られるプロファイルが滑らかに結合されている様子がわかります。

まとめ

画像処理において強度プロファイルを作成する方法について、scikit-imageライブラリのprofile_line()関数を使用して解説しました。この関数を用いることで、任意の直線上のプロファイルを簡単に取得できます。また、この関数の主要な引数について説明し、実装例を通じて使い方を示しました。

プロファイルの作成や解析は画像の性質を理解する重要な手法ですので、本記事で紹介した方法が、効率的かつ正確な画像解析を行うための一助となれば幸いです。

参考資料