7

ここに画像の説明を入力

オペレーターは、各ピークの位置と幅を把握してスペクトルを調べ、スペクトルが属する部分を判断していました。新しい方法では、画像はカメラによってスクリーンに取り込まれます。また、各バンドの幅はプログラムで計算する必要があります。

古いシステム: 分光器 -> 人間の目 新しいシステム: 分光器 -> カメラ -> プログラム

おおよそのX軸位置が与えられた場合、各バンドの幅を計算するための良い方法は何ですか? このタスクは以前は目で完全に実行されていましたが、現在はプログラムで実行する必要がありますか?

詳細が不足している場合は申し訳ありませんが、不足しています。

前のグラフを生成したプログラムのリスト。私はそれが関連していることを願っています:

import Image
from scipy import *
from scipy.optimize import leastsq

# Load the picture with PIL, process if needed
pic         = asarray(Image.open("spectrum.jpg"))

# Average the pixel values along vertical axis
pic_avg     = pic.mean(axis=2)
projection  = pic_avg.sum(axis=0)

# Set the min value to zero for a nice fit
projection /= projection.mean()
projection -= projection.min()

#print projection

# Fit function, two gaussians, adjust as needed
def fitfunc(p,x):
    return p[0]*exp(-(x-p[1])**2/(2.0*p[2]**2)) + \
        p[3]*exp(-(x-p[4])**2/(2.0*p[5]**2))
errfunc = lambda p, x, y: fitfunc(p,x)-y

# Use scipy to fit, p0 is inital guess
p0 = array([0,20,1,0,75,10])
X  = xrange(len(projection))
p1, success = leastsq(errfunc, p0, args=(X,projection))
Y = fitfunc(p1,X)

# Output the result
print "Mean values at: ", p1[1], p1[4]

# Plot the result
from pylab import *
#subplot(211)
#imshow(pic)
#subplot(223)
#plot(projection)
#subplot(224)
#plot(X,Y,'r',lw=5)
#show()

subplot(311)
imshow(pic)
subplot(312)
plot(projection)
subplot(313)
plot(X,Y,'r',lw=5)
show()
4

3 に答える 3

3

おおよその開始点が与えられると、この点に最も近い極大値を見つける単純なアルゴリズムを使用できます。あなたの適切なコードはすでにそれを行っている可能性があります（あなたがそれをうまく使っているかどうかはわかりませんでした）。

ユーザーが指定した開始点からの単純なピーク検出を示すコードを次に示します。

#!/usr/bin/env python
from __future__ import division
import numpy as np
from matplotlib import pyplot as plt

# Sample data with two peaks: small one at t=0.4, large one at t=0.8
ts = np.arange(0, 1, 0.01)
xs = np.exp(-((ts-0.4)/0.1)**2) + 2*np.exp(-((ts-0.8)/0.1)**2)

# Say we have an approximate starting point of 0.35
start_point = 0.35

# Nearest index in "ts" to this starting point is...
start_index = np.argmin(np.abs(ts - start_point))

# Find the local maxima in our data by looking for a sign change in
# the first difference
# From http://stackoverflow.com/a/9667121/188535
maxes = (np.diff(np.sign(np.diff(xs))) < 0).nonzero()[0] + 1

# Find which of these peaks is closest to our starting point
index_of_peak = maxes[np.argmin(np.abs(maxes - start_index))]

print "Peak centre at: %.3f" % ts[index_of_peak]

# Quick plot showing the results: blue line is data, green dot is
# starting point, red dot is peak location
plt.plot(ts, xs, '-b')
plt.plot(ts[start_index], xs[start_index], 'og')
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or')
plt.show()

この方法は、ピークへの上昇が開始点から完全にスムーズである場合にのみ機能します。これがノイズに対する耐性を高める必要がある場合、私は使用していませんが、PyDSToolが役立つようです。このSciPy の投稿では、ノイズの多いデータ セットで 1D ピークを検出するために SciPy を使用する方法について詳しく説明しています。

この時点で、ピークの中心を見つけたとします。幅については、使用できる方法がいくつかありますが、最も簡単なのはおそらく「半値全幅」(FWHM) です。繰り返しますが、これは単純であるため脆弱です。2 つのピークが近い場合や、ノイズの多いデータの場合は壊れます。

FWHM はまさにその名前が示すとおりです。ピークの幅が最大値の中間にあることを示します。これを行うコードを次に示します (上記の続きです)。

# FWHM...
half_max = xs[index_of_peak]/2

# This finds where in the data we cross over the halfway point to our peak. Note
# that this is global, so we need an extra step to refine these results to find
# the closest crossovers to our peak.

# Same sign-change-in-first-diff technique as above
hm_left_indices = (np.diff(np.sign(np.diff(np.abs(xs[:index_of_peak] - half_max)))) > 0).nonzero()[0] + 1
# Add "index_of_peak" to result because we cut off the left side of the data!
hm_right_indices = (np.diff(np.sign(np.diff(np.abs(xs[index_of_peak:] - half_max)))) > 0).nonzero()[0] + 1 + index_of_peak

# Find closest half-max index to peak
hm_left_index = hm_left_indices[np.argmin(np.abs(hm_left_indices - index_of_peak))]
hm_right_index = hm_right_indices[np.argmin(np.abs(hm_right_indices - index_of_peak))]

# And the width is...    
fwhm = ts[hm_right_index] - ts[hm_left_index]

print "Width: %.3f" % fwhm

# Plot to illustrate FWHM: blue line is data, red circle is peak, red line
# shows FWHM
plt.plot(ts, xs, '-b')
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or')
plt.plot(
    [ts[hm_left_index], ts[hm_right_index]],
    [xs[hm_left_index], xs[hm_right_index]], '-r')
plt.show()

半値全幅である必要はありません— あるコメンターが指摘しているように、オペレータのピーク検出の通常のしきい値がどこにあるかを把握し、それをプロセスのこのステップのアルゴリズムに変えることができます。

より堅牢な方法は、ガウス曲線 (または独自のモデル) をピークを中心とするデータのサブセットに適合させることです。たとえば、一方の極小値から他方の極小値まで、幅を計算するためのその曲線のパラメーター (例: シグマ)。

これが大量のコードであることは承知していますが、インデックス検索関数を除外して「自分の作業を表示」することは意図的に避けています。

これが少なくとも、特定のセットにより適したものを考え出すための良い出発点になることを願っています.

于 2012-05-26T09:09:36.013 に答える
2

パーティーに遅れましたが、将来この質問に出くわす人のために...

眼球運動データはこれと非常によく似ています。私は、Nystrom + Holmqvist、2010で使用されたアプローチに基づいています。Savitsky-Golay フィルター ( scipy.signal.savgol_filterscipy v0.14+) を使用してデータを平滑化して、大きなピークをそのまま維持しながら低レベルのノイズの一部を取り除きます。著者は、2 の次数と約 2 倍のウィンドウ サイズを使用することを推奨しています検出できるようにしたい最小ピークの幅。特定の y 値を超えるすべての値を任意に削除することで、バンドがどこにあるかを見つけることができます (それらをnumpy.nan）。次に、残りの（nan）平均と（nan）標準偏差を取り、平均+ [パラメータ] * std（論文では6を使用していると思います）より大きいすべての値を削除します。データポイントを削除しなくなるまで繰り返しますが、データによっては、[パラメーター] の特定の値が安定しない場合があります。次にnumpy.isnan()、イベントと非イベントnumpy.diff()を検索し、各イベントの開始と終了 (それぞれ -1 と 1 の値) を見つけるために使用します。さらに正確な開始点と終了点を取得するには、データに沿って各開始点から逆方向にスキャンし、各終了点から順方向にスキャンして、平均 + [別のパラメーター]*std よりも小さい値を持つ最も近いローカル最小値を見つけることができます (3 を使用すると思います)。紙で）。次に、各開始点と終了点の間のデータ ポイントを数えるだけです。

これは、その二重ピークでは機能しません。そのためには、いくつかの外挿を行う必要があります。

于 2014-11-14T00:46:39.150 に答える
0

最良の方法は、一連の方法を人間の結果と統計的に比較することかもしれません。

多種多様なデータと多種多様な測定推定値（さまざまなしきい値での幅、さまざまなしきい値を超える領域、さまざまなしきい値の選択方法、2次モーメント、さまざまな次数の多項式カーブフィット、パターンマッチングなど）を取得し、これらを比較します。同じデータセットの人間の測定値に対する推定。専門家の人間の結果と最もよく相関する推定方法を選択してください。または、さまざまな高さのそれぞれに最適な方法、他のピークからのさまざまな分離など、いくつかの方法を選択することもできます。

于 2012-05-26T21:19:57.027 に答える