这些光谱带过去是用眼睛来判断的,如何以编程方式进行？

Question

操作员用于检查光谱,知道每个峰的位置和宽度,并判断光谱所属的部分.以新的方式,图像由相机捕获到屏幕.并且必须以编程方式计算每个波段的宽度.

旧系统:分光镜 - >人眼新系统:分光镜 - >摄像头 - >程序

计算每个波段宽度的好方法是什么,给出它们近似的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()

Answer 1

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

\n\n

下面是一些代码，演示了从用户给定的起点进行简单的峰值查找：

\n\n

#!/usr/bin/env python\nfrom __future__ import division\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\n# Sample data with two peaks: small one at t=0.4, large one at t=0.8\nts = np.arange(0, 1, 0.01)\nxs = np.exp(-((ts-0.4)/0.1)**2) + 2*np.exp(-((ts-0.8)/0.1)**2)\n\n# Say we have an approximate starting point of 0.35\nstart_point = 0.35\n\n# Nearest index in "ts" to this starting point is...\nstart_index = np.argmin(np.abs(ts - start_point))\n\n# Find the local maxima in our data by looking for a sign change in\n# the first difference\n# From http://stackoverflow.com/a/9667121/188535\nmaxes = (np.diff(np.sign(np.diff(xs))) < 0).nonzero()[0] + 1\n\n# Find which of these peaks is closest to our starting point\nindex_of_peak = maxes[np.argmin(np.abs(maxes - start_index))]\n\nprint "Peak centre at: %.3f" % ts[index_of_peak]\n\n# Quick plot showing the results: blue line is data, green dot is\n# starting point, red dot is peak location\nplt.plot(ts, xs, \'-b\')\nplt.plot(ts[start_index], xs[start_index], \'og\')\nplt.plot(ts[index_of_peak], xs[index_of_peak], \'or\')\nplt.show()\n

\n\n

仅当从起点开始完全顺利地攀登山顶时，此方法才有效。如果这需要对噪音更有弹性，我还没有使用它，但PyDSTool似乎可能会有所帮助。这篇SciPy 帖子详细介绍了如何使用它来检测噪声数据集中的一维峰值。

\n\n

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

\n\n

FWHM 正是顾名思义：您可以找到峰在最大值一半处的宽度。这是一些执行此操作的代码（它只是从上面继续）：

\n\n

# FWHM...\nhalf_max = xs[index_of_peak]/2\n\n# This finds where in the data we cross over the halfway point to our peak. Note\n# that this is global, so we need an extra step to refine these results to find\n# the closest crossovers to our peak.\n\n# Same sign-change-in-first-diff technique as above\nhm_left_indices = (np.diff(np.sign(np.diff(np.abs(xs[:index_of_peak] - half_max)))) > 0).nonzero()[0] + 1\n# Add "index_of_peak" to result because we cut off the left side of the data!\nhm_right_indices = (np.diff(np.sign(np.diff(np.abs(xs[index_of_peak:] - half_max)))) > 0).nonzero()[0] + 1 + index_of_peak\n\n# Find closest half-max index to peak\nhm_left_index = hm_left_indices[np.argmin(np.abs(hm_left_indices - index_of_peak))]\nhm_right_index = hm_right_indices[np.argmin(np.abs(hm_right_indices - index_of_peak))]\n\n# And the width is...    \nfwhm = ts[hm_right_index] - ts[hm_left_index]\n\nprint "Width: %.3f" % fwhm\n\n# Plot to illustrate FWHM: blue line is data, red circle is peak, red line\n# shows FWHM\nplt.plot(ts, xs, \'-b\')\nplt.plot(ts[index_of_peak], xs[index_of_peak], \'or\')\nplt.plot(\n    [ts[hm_left_index], ts[hm_right_index]],\n    [xs[hm_left_index], xs[hm_right_index]], \'-r\')\nplt.show()\n

\n\n

它不必是半宽最大值 \xe2\x80\x94 的全宽，您可以尝试找出操作员峰值检测的正常阈值在哪里，并将其转化为一种算法对于该过程的这一步。

\n\n

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

\n\n

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

\n\n

希望这至少能为您提供一个良好的起点，以提出更适合您的特定集合的东西。

\n