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

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

# 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|ピーク検出の解像度 <br> 値が大きいほど荒くピークを探索|
|order|ピーク検出の閾値 <br> 最大ピーク値のorder倍以下のピークを排除|

|返り値 |説明  |
|---|---|
|peaks_index|ピークが格納されている要素のインデックス|

- シミュレーション  
以下の2次元配列データからピークを検出してみます.
右図は左図を上から見たものです.
<img width="936" alt="スクリーンショット 2019-04-04 19.29.18.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/c331fb0d-5601-937e-7cd5-4adab6346adc.png">

ピーク検出の結果
<img width="1007" alt="スクリーンショット 2019-04-04 19.23.38.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/b898f18d-8c47-9cc3-1b2d-7c3aa7c60af4.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. 低値ピークの排除

今回は、以下のようなデータを例に説明します.
元配列
<img width="312" alt="スクリーンショット 2019-04-04 18.15.44.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/e92bdcca-11b9-8946-23f4-b158efb654c4.png">
上の配列を可視化したもの
<img width="939" alt="スクリーンショット 2019-04-04 19.27.08.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/6c0e7607-0992-494f-c22b-7fc8936c7d20.png">


1. 最大値フィルタによる処理
まず、`filter_size`で設定した大きさの最大値フィルタでデータをのっぺりとさせます.
3×3サイズの最大値フィルタで処理後
<img width="307" alt="スクリーンショット 2019-04-04 18.17.56.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/d2981164-9cd9-7c11-6794-bb3dddb5c7d9.png">
<img width="928" alt="スクリーンショット 2019-04-04 19.28.34.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/bf6dd1e9-644b-93ce-0b65-0e93dddd024b.png">


2. マスク処理
元の配列と同じ値を持つ要素だけ残し、それ以外の要素を０にします.
そうすることで、周りより値の高い要素だけが残ります.
<img width="313" alt="スクリーンショット 2019-04-04 18.32.18.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/9caa1969-4479-7fe3-3067-6fac2ce53fa8.png">
<img width="943" alt="スクリーンショット 2019-04-04 19.25.33.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/04feab42-61b6-e7dc-c6cf-a7cf59115bd9.png">


3. 低値ピークの排除
`order`で設定した値に応じて、小さい値のピークを排除します.
<img width="304" alt="スクリーンショット 2019-04-04 18.56.16.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/a37185ee-d039-4265-c0ac-6c15c2725e7a.png">
<img width="935" alt="スクリーンショット 2019-04-04 19.26.25.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/e1876366-390e-2929-379c-1e0116f92a48.png">

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

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

# おまけ
## １次元配列データのピーク検出
[scipy.signal.argrelmax](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.argrelmax.html)を使用すればできます.
引数の`order`の値によって検出するピークの閾値を決定します.
値が低いほど、敏感にピークを検出します.

- シミュレーション
<img width="986" alt="スクリーンショット 2019-04-04 16.08.35.png" src="https://qiita-image-store.s3.amazonaws.com/0/301489/c00052c1-0852-d292-5d89-581deeb2349b.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()
```