这些频段过去是用肉眼来判断的，如何通过程序来判断呢？答案

【问题标题】：These spectrum bands used to be judged by eye, how to do it programmatically?这些频段过去是用肉眼来判断的，如何通过程序来判断呢？
【发布时间】：2012-06-01 15:15:10
【问题描述】：

操作员用于检查光谱，知道每个峰的位置和宽度，并判断光谱所属的片段。在新的方式中，图像由相机捕捉到屏幕上。并且必须以编程方式计算每个波段的宽度。

旧系统：分光镜 -> 人眼新系统：分光镜 -> 相机 -> 程序

根据 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()

【问题讨论】：

找到人眼无法区分颜色和背景的阈值。它应该是相当恒定的。然后，只需对您的数据点设置阈值，使人眼“看到”高于该阈值的数据点，然后对这些数据点进行聚类并找到每个聚类的宽度。
也许对曲线求导很容易找到峰值和斜率？条带的宽度不是很恒定吗？据我了解，振幅和位置会有所不同。
@Blender，你能详细说明你的观点吗？
本底噪声是否相当平坦？您如何（通过肉眼）确定什么是高峰，什么是太低？例如，大约 660-670 处的那个光点是双峰，还是不会被计算在内？ 740 的时候呢？
@detly，我并不完全了解，但峰的位置会事先知道。是的，如果两个已知的山峰彼此太近，我会担心。

标签： python image-processing signal-processing

【解决方案1】：

给定一个近似的起点，您可以使用一个简单的算法来找到最接近该点的局部最大值。您的拟合代码可能已经这样做了（我不确定您是否成功使用它）。

以下代码演示了从用户给定的起点进行简单的峰值查找：

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

因此假设此时您已找到峰的中心。现在对于宽度：您可以使用几种方法，但最简单的可能是“半高全宽”（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()

它不必是最大一半的全宽 — 正如一位评论者指出的那样，您可以尝试找出操作员的峰值检测正常阈值在哪里，然后将其转换为进入该过程的这一步的算法。

一种更稳健的方法可能是将高斯曲线（或您自己的模型）拟合到以峰值为中心的数据子集（例如，从一侧的局部最小值到另一侧的局部最小值）并使用该曲线的参数之一（例如 sigma）来计算宽度。

我意识到这是很多代码，但我故意避免将索引查找函数分解为更多“展示我的工作”，当然绘图函数只是为了演示。

希望这至少能给你一个很好的起点，让你想出更适合你的特定系列的东西。

【讨论】：

如果“基础”低于阈值（在这种情况下为峰值的 0.5），您能否添加一个引发异常的修改？
@aitchnyu 好吧，我在这里做了一点过度简化：我刚刚说的是总高度的 0.5，但 FWHM 最好测量为峰高度的一半 在底部之上（或本底噪声）。所以基数低于峰值的 0.5 是没有意义的。
是的，如果算法遇到不正确的峰值，我希望算法抛出异常：峰值不会单调下降到 0.5 个峰值（请参见我图中的双峰）。我只需要宽度，幅度会因仪器的特性而改变。 “本底噪声”在理想情况下应该是无声的，对应于黑色的强度。
@aitchnyu 啊，我明白了。我没有时间自己做，但是您可以实现一些查看峰值范围内的第一个差异的东西。它应该在峰的左侧严格为正，在右侧严格为负。根据定义，它在峰值处为零，但也许您也应该允许任一侧的值为零。
所以第一个差异的hm_left_index:index_of_peak 必须全为正，index_of_peak:hm_right_index 必须全为负，对吧？

【解决方案2】：

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

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

【讨论】：

【解决方案3】：

聚会迟到了，但是对于以后遇到这个问题的任何人...

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

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

【讨论】：