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 文章详细介绍了如何使用它来检测嘈杂数据集中的一维峰值。

所以假设此时你已经找到了峰值的中心。现在对于宽度:您可以使用多种方法,但最简单的可能是“半高全宽”(FWHM)。同样,这很简单,因此很脆弱。它会因接近双峰或噪声数据而中断。

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_filter在 scipy v0.14+ 中)对数据进行平滑处理,以消除一些低级噪声,同时保持大峰完好无损——作者建议使用 2 阶和大约两倍的窗口大小您希望能够检测到的最小峰的宽度。您可以通过任意删除高于某个 y 值的所有值来找到波段的位置(将它们设置为numpy.nan)。然后取余数的 (nan)mean 和 (nan) 标准差,并删除所有大于均值 + [parameter]*std 的值(我认为他们在论文中使用 6)。迭代直到您不删除任何数据点 - 但根据您的数据,[parameter] 的某些值可能不稳定。然后用于numpy.isnan()查找事件与非事件,并numpy.diff()查找每个事件的开始和结束(值分别为 -1 和 1)。为了获得更准确的起点和终点,您可以从每个起点向后扫描数据并从每一端向前扫描,以找到最近的局部最小值,其值小于 mean + [another parameter]*std(我认为他们使用 3在论文中)。然后你只需要计算每个开始和结束之间的数据点。

这不适用于双峰;你必须为此做一些推断。

于 2014-11-14T00:46:39.150 回答
0

最好的方法可能是将一堆方法与人类结果进行统计比较。

您将获取大量数据和大量测量估计(各种阈值处的宽度、各种阈值以上的区域、不同的阈值选择方法、二阶矩、各种度数的多项式曲线拟合、模式匹配等)并比较这些对同一数据集的人类测量值的估计。选择与专家人类结果最相关的估计方法。或者可以选择几种方法,对于各种高度中的每一种,对于与其他峰的各种分离,等等都是最好的。

于 2012-05-26T21:19:57.027 回答