46
41

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

Pythonでピーク検出(1次元配列データと2次元配列データ)

Last updated at Posted at 2019-04-04

はじめに

Pythonで2次元配列データのピーク検出のライブラリが見当たらなかったので、自分で関数を作ってみました.

2次元配列データのピーク検出

結論から言うと、以下の関数を使用すればピークを検出できます.
関数が何をやっているかについては2次元配列データのピーク検出の仕組みで説明しています.

# ピーク検出関数
def detect_peaks(image, filter_size=3, order=0.5):
    local_max = maximum_filter(image, footprint=np.ones((filter_size, filter_size)), mode='constant')
    detected_peaks = np.ma.array(image, mask=~(image == local_max))

    # 小さいピーク値を排除(最大ピーク値のorder倍以下のピークは排除)
    temp = np.ma.array(detected_peaks, mask=~(detected_peaks >= detected_peaks.max() * order))
    peaks_index = np.where((temp.mask != True))
    return peaks_index
引数 説明
image 入力データ(2D numpy array)
filter_size ピーク検出の解像度
値が大きいほど荒くピークを探索
order ピーク検出の閾値
最大ピーク値のorder倍以下のピークを排除
返り値 説明
peaks_index ピークが格納されている要素のインデックス
  • シミュレーション
    以下の2次元配列データからピークを検出してみます.
    右図は左図を上から見たものです.
スクリーンショット 2019-04-04 19.29.18.png

ピーク検出の結果
スクリーンショット 2019-04-04 19.23.38.png

  • サンプルコード
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import maximum_filter


# ピーク検出関数
def detect_peaks(image, filter_size=3, order=0.5):
    local_max = maximum_filter(image, footprint=np.ones((filter_size, filter_size)), mode='constant')
    detected_peaks = np.ma.array(image, mask=~(image == local_max))

    # 小さいピーク値を排除(最大ピーク値のorder倍のピークは排除)
    temp = np.ma.array(detected_peaks, mask=~(detected_peaks >= detected_peaks.max() * order))
    peaks_index = np.where((temp.mask != True))
    return peaks_index

# データ読み込み(好きな2d arrayを読み込んでください)
Map_val = np.loadtxt('Map2_val2.txt')
Map_val = np.rot90(Map_val, k=-2)

# ---------- ピーク検出 ---------- #
maxid1 = detect_peaks(Map_val, order=0.3)
maxid2 = detect_peaks(Map_val, order=0.5)


# ---------- 描画 ---------- #
Fig = plt.figure(figsize=(12, 6))
Map1 = Fig.add_subplot(121)
Map2 = Fig.add_subplot(122)

Map1.set_title('order=0.3')
Map1.imshow(Map_val.T, interpolation='spline36', cmap="rainbow")
for i in range(len(maxid1[0])):
    Map1.scatter(maxid1[0][i], maxid1[1][i], color='black')
    Map1.text(maxid1[0][i], maxid1[1][i], 'PEAK!!!')

Map2.set_title('order=0.5')
Map2.imshow(Map_val.T, interpolation='spline36', cmap="rainbow")
for i in range(len(maxid2[0])):
    Map2.scatter(maxid2[0][i], maxid2[1][i], color='black')
    Map2.text(maxid2[0][i], maxid2[1][i], 'PEAK!!!')

plt.show()

2次元配列データのピーク検出の仕組み

簡単に言うと周りより値が高いところをピークとして検出しています.
その後、小さい値のピークを排除しているだけです.
「周り」をどれぐらいまで考慮するかのパラメータをfilter_sizeで、「小さい」っていう尺度をorderで決定しています.
具体的には、以下の処理によりピークを検出しています.

  1. 最大値フィルタによる処理
  2. マスク処理
  3. 低値ピークの排除

今回は、以下のようなデータを例に説明します.
元配列
スクリーンショット 2019-04-04 18.15.44.png
上の配列を可視化したもの
スクリーンショット 2019-04-04 19.27.08.png

  1. 最大値フィルタによる処理
    まず、filter_sizeで設定した大きさの最大値フィルタでデータをのっぺりとさせます.
    3×3サイズの最大値フィルタで処理後
スクリーンショット 2019-04-04 18.17.56.png スクリーンショット 2019-04-04 19.28.34.png
  1. マスク処理
    元の配列と同じ値を持つ要素だけ残し、それ以外の要素を0にします.
    そうすることで、周りより値の高い要素だけが残ります.
スクリーンショット 2019-04-04 18.32.18.png スクリーンショット 2019-04-04 19.25.33.png
  1. 低値ピークの排除
    orderで設定した値に応じて、小さい値のピークを排除します.
スクリーンショット 2019-04-04 18.56.16.png スクリーンショット 2019-04-04 19.26.25.png

これでピークを検出することができました.

最後に

 今回は「値が周りより大きい所」and「値が全体の平均より大幅に大きい所」をピークとして検出する関数を作りました.
しかし、「値が周りより大きい所」and「値が周りの平均値より大幅に大きい所(周りとの変化が大きい所)」をピークと考えることもでき、その場合は「3.低値ピークの排除」のところの処理が変わりますので、興味ある方はぜひ実装して見てください(平均値フィルタを使えば実装できると思います).
 またピーク検出について調べていると、そもそもピークをどう定義するかやノイズが大きい場合のノイズとピークの違いなど、なかなか奥が深いことがわかりました.

おまけ

1次元配列データのピーク検出

scipy.signal.argrelmaxを使用すればできます.
引数のorderの値によって検出するピークの閾値を決定します.
値が低いほど、敏感にピークを検出します.

  • シミュレーション
スクリーンショット 2019-04-04 16.08.35.png
  • サンプルコード
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import signal

# 人工データの作成
N = 100
x = np.linspace(-1, 1, N)
theta = 2 * math.pi * ((x + 1.0) / 2.0)
y = np.sin(theta) + np.random.randn(N)*0.05

# ピーク値のインデックスを取得
maxid = signal.argrelmax(y, order=5) # 最大値

# 描画
Fig = plt.figure(figsize=(6, 6))
Map1 = Fig.add_subplot(111)
Map1.plot(x, y)
Map1.scatter(x[maxid[0]], y[maxid[0]])
for i in range(len(maxid[0])):
    Map1.text(x[maxid[0][i]], y[maxid[0][i]], 'PEAK!!!')

plt.show()
46
41
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
46
41

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?