2D阵列中的峰值检测

我正在帮助一家兽医诊所测量狗爪下的压力。我使用Python进行数据分析，现在我不得不试图将爪子分成（解剖学）子区域。

我制作了每个爪子的2D阵列，它由爪子随时间加载的每个传感器的最大值组成。这是一个爪子的例子，我用Excel绘制了我想要“检测”的区域。这些是传感器周围的2×2个盒子，具有局部最大值，它们一起具有最大的总和。

所以我尝试了一些实验并决定只查找每列和每行的最大值（由于爪子的形状，不能在一个方向上查看）。这似乎可以很好地“检测”单独脚趾的位置，但它也标记了相邻的传感器。

那么告诉Python哪些最大值是我想要的最好的方法是什么？

注意：2x2正方形不能重叠，因为它们必须是单独的脚趾！

我也采用2x2作为方便，欢迎任何更高级的解决方案，但我只是一个人类运动科学家，所以我既不是真正的程序员也不是数学家，所以请保持“简单”。

这是一个可以用np.loadtxt加载的版本

结果

所以我尝试了@jextee的解决方案（见下面的结果）。正如你所看到的，它在前爪上很有效，但后腿的效果不太好。

更具体地说，它无法识别出第四个脚趾的小峰值。这显然是循环看起来自上而下朝向最低值的事实所固有的，而不考虑这是什么。

有谁知道如何调整@jextee的算法，以便它也可以找到第4个脚趾？

由于我还没有处理任何其他试验，我不能提供任何其他样品。但我之前提供的数据是每只爪子的平均值。该文件是一个数组，其最大数据为9个爪子，它们与盘子接触的顺序。

该图像显示了它们如何在空间上展开。

更新：

我已经为任何感兴趣的人建立了一个博客 ，我已经设置了一个包含所有原始测量值的SkyDrive。所以对于要求更多数据的人来说：给你更大的力量！

新更新：

所以在我得到关于爪子检测和爪子分类的问题的帮助后，我终于能够检查每个爪子的脚趾检测！事实证明，除了像我自己的例子中那样大小的爪子之外，它在任何东西上都不能很好地工作。事后看来，任意选择2x2是我自己的错。

这是一个错误的例子：钉子被识别为脚趾，“脚跟”如此宽，它被识别两次！

爪子太大，因此在没有重叠的情况下采用2x2尺寸会导致一些脚趾被检测到两次。相反，在小型犬中，它经常无法找到第五个脚趾，我怀疑它是由2x2区域太大引起的。

在对我的所有测量结果进行了当前的解决方案之后，我得出了令人吃惊的结论：几乎所有的小型犬都没有找到第5个脚趾，并且对于大型犬的50％以上的影响它会发现更多！

显然我需要改变它。我自己的猜测是，对于小型犬而言，将neighborhood的大小改为小型犬，对于大型犬则更大。但是generate_binary_structure不会让我改变数组的大小。

因此，我希望其他人有更好的建议来定位脚趾，也许脚趾区域尺寸与爪子尺寸一致？

#1楼

我相信你现在已经足够了，但我不禁建议使用k-means聚类方法。 k-means是一种无监督的聚类算法，可以获取数据（在任意数量的维度中 - 我碰巧在3D中这样做）并将其排列成具有不同边界的k个聚类。这里很好，因为你确切地知道这些犬（应该）有多少脚趾。

另外，它在Scipy中实现，这非常好（ http://docs.scipy.org/doc/scipy/reference/cluster.vq.html ）。

以下是它可以在空间上解析3D集群的示例：

你想要做的是有点不同（2D并包括压力值），但我仍然认为你可以试一试。

#2楼

我只想了几个想法：

采取扫描的梯度（衍生），看看是否消除了错误的调用
取最大的局部最大值

您可能还想看看OpenCV ，它有一个相当不错的Python API，并且可能有一些您觉得有用的功能。

#3楼

也许一个天真的方法就足够了：在你的飞机上建立一个所有2x2方格的列表，按它们的总和（按降序排列）排序。

首先，在“爪子列表”中选择值最高的方格。然后，迭代地选择4个与之前找到的任何方格都不相交的下一个最佳方格。

#4楼

这是一个想法：你计算图像的（离散）拉普拉斯算子。我希望它在最大值时是（负的）大，其方式比原始图像更具戏剧性。因此，最大值可能更容易找到。

这是另一个想法：如果你知道高压点的典型大小，你可以先用相同大小的高斯卷积它来平滑你的图像。这可能会让您处理更简单的图像。

#5楼

如果您一步一步地执行操作：首先找到全局最大值，如果需要，处理给定值的周围点，然后将找到的区域设置为零，并重复下一个。

#6楼

解

数据文件： paw.txt 。源代码：

from scipy import *
from operator import itemgettern = 5  # how many fingers are we looking ford = loadtxt("paw.txt")
width, height = d.shape# Create an array where every element is a sum of 2x2 squares.fourSums = d[:-1,:-1] + d[1:,:-1] + d[1:,1:] + d[:-1,1:]# Find positions of the fingers.# Pair each sum with its position number (from 0 to width*height-1),pairs = zip(arange(width*height), fourSums.flatten())# Sort by descending sum value, filter overlapping squaresdef drop_overlapping(pairs):no_overlaps = []def does_not_overlap(p1, p2):i1, i2 = p1[0], p2[0]r1, col1 = i1 / (width-1), i1 % (width-1)r2, col2 = i2 / (width-1), i2 % (width-1)return (max(abs(r1-r2),abs(col1-col2)) >= 2)for p in pairs:if all(map(lambda prev: does_not_overlap(p,prev), no_overlaps)):no_overlaps.append(p)return no_overlapspairs2 = drop_overlapping(sorted(pairs, key=itemgetter(1), reverse=True))# Take the first n with the heighest valuespositions = pairs2[:n]# Print resultsprint d, "\n"for i, val in positions:row = i / (width-1)column = i % (width-1)print "sum = %f @ %d,%d (%d)" % (val, row, column, i)print d[row:row+2,column:column+2], "\n"

输出没有重叠方块。似乎在您的示例中选择了相同的区域。

一些评论

棘手的部分是计算所有2x2平方的总和。我假设你需要所有这些，所以可能会有一些重叠。我使用切片从原始2D数组中剪切第一个/最后一个列和行，然后将它们重叠在一起并计算总和。

为了更好地理解它，对3x3阵列进行成像：

>>> a = arange(9).reshape(3,3) ; a
array([[0, 1, 2],[3, 4, 5],[6, 7, 8]])

然后你可以拿它的切片：

>>> a[:-1,:-1]
array([[0, 1],[3, 4]])
>>> a[1:,:-1]
array([[3, 4],[6, 7]])
>>> a[:-1,1:]
array([[1, 2],[4, 5]])
>>> a[1:,1:]
array([[4, 5],[7, 8]])

现在想象一下，将它们叠加在另一个之上，并将元素加在相同位置。这些总和与2x2正方形的总和相同，左上角位于同一位置：

>>> sums = a[:-1,:-1] + a[1:,:-1] + a[:-1,1:] + a[1:,1:]; sums
array([[ 8, 12],[20, 24]])

如果总和超过2x2平方，则可以使用max查找最大值，或sort或sorted以查找峰值。

为了记住峰值的位置，我将每个值（总和）与其在序列阵列中的序数位置耦合（参见zip ）。然后我打印结果时再次计算行/列位置。

笔记

我允许2x2方格重叠。编辑版本过滤掉其中一些版本，使得结果中仅出现非重叠的正方形。

选择手指（一个想法）

另一个问题是如何从所有峰中选择可能是手指的东西。我有一个想法可能会或可能不会奏效。我现在没时间实现它，所以只是伪代码。

我注意到如果前手指几乎保持完美的圆形，则后指应位于该圆圈的内侧。而且，前指或多或少地等间隔。我们可能会尝试使用这些启发式属性来检测手指。

伪代码：

select the top N finger candidates (not too many, 10 or 12)
consider all possible combinations of 5 out of N (use itertools.combinations)
for each combination of 5 fingers:for each finger out of 5:fit the best circle to the remaining 4=> position of the center, radiuscheck if the selected finger is inside of the circlecheck if the remaining four are evenly spread(for example, consider angles from the center of the circle)assign some cost (penalty) to this selection of 4 peaks + a rear finger(consider, probably weighted:circle fitting error,if the rear finger is inside,variance in the spreading of the front fingers,total intensity of 5 peaks)
choose a combination of 4 peaks + a rear peak with the lowest penalty

这是一种蛮力的方法。如果N相对较小，那么我认为它是可行的。对于N = 12，有C_12 ^ 5 = 792种组合，有5种方式可以选择一个后指，所以3960个案例来评估每个爪子。

#7楼

好吧，这里有一些简单而且效率不高的代码，但对于这个数据集的大小，这很好。

import numpy as np
grid = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0.4,0.4,0.4,0,0,0],[0,0,0,0,0.4,1.4,1.4,1.8,0.7,0,0,0,0,0],[0,0,0,0,0.4,1.4,4,5.4,2.2,0.4,0,0,0,0],[0,0,0.7,1.1,0.4,1.1,3.2,3.6,1.1,0,0,0,0,0],[0,0.4,2.9,3.6,1.1,0.4,0.7,0.7,0.4,0.4,0,0,0,0],[0,0.4,2.5,3.2,1.8,0.7,0.4,0.4,0.4,1.4,0.7,0,0,0],[0,0,0.7,3.6,5.8,2.9,1.4,2.2,1.4,1.8,1.1,0,0,0],[0,0,1.1,5,6.8,3.2,4,6.1,1.8,0.4,0.4,0,0,0],[0,0,0.4,1.1,1.8,1.8,4.3,3.2,0.7,0,0,0,0,0],[0,0,0,0,0,0.4,0.7,0.4,0,0,0,0,0,0]])arr = []
for i in xrange(grid.shape[0] - 1):for j in xrange(grid.shape[1] - 1):tot = grid[i][j] + grid[i+1][j] + grid[i][j+1] + grid[i+1][j+1]arr.append([(i,j),tot])best = []arr.sort(key = lambda x: x[1])for i in xrange(5):best.append(arr.pop())badpos = set([(best[-1][0][0]+x,best[-1][0][1]+y)for x in [-1,0,1] for y in [-1,0,1] if x != 0 or y != 0])for j in xrange(len(arr)-1,-1,-1):if arr[j][0] in badpos:arr.pop(j)for item in best:print grid[item[0][0]:item[0][0]+2,item[0][1]:item[0][1]+2]

我基本上只是制作一个数组，其中左上角的位置和每个2x2平方的总和，并按总和进行排序。然后我拿出2x2平方的最高和的争用，把它放在best阵列中，并删除所有其他2x2正方形，使用这个刚删除的2x2平方的任何部分。

它似乎工作正常，除了最后一个爪子（你的第一张照片中最右边的一个最小的一个），事实证明，还有另外两个符合条件的2x2正方形和一个更大的总和（他们有一个相等的总和彼此）。其中一个仍然从你的2x2方格中选择一个方格，但另一个方向偏离左边。幸运的是，幸运的是，我们看到选择更多你想要的，但这可能需要一些其他的想法来获得你真正想要的东西。

#8楼

这是图像注册问题。总体战略是：

有一个已知的例子，或者某种先前的数据。
使您的数据适合示例，或使示例适合您的数据。
如果您的数据首先大致一致，这会有所帮助。

这是一个粗略而准备好的方法 ，“可能有用的最蠢的事情”：

从大约您期望的位置开始五个脚趾坐标。
每一个，迭代地爬到山顶。即，给定当前位置，如果其值大于当前像素，则移动到最大相邻像素。当你的脚趾坐标停止移动时停止。

要抵消方向问题，您可以为基本方向（北，东北等）设置8个左右的初始设置。单独运行每一个并丢弃任何结果，其中两个或多个脚趾最终在同一像素。我会更多地考虑这个问题，但是这种事情仍然在图像处理中进行研究 - 没有正确的答案！

稍微复杂的想法:(加权）K均值聚类。 没那么糟糕。

从五个脚趾坐标开始，但现在这些是“集群中心”。

然后迭代直到收敛：

将每个像素分配给最近的群集（只为每个群集创建一个列表）。
计算每个群集的质心。对于每个聚类，它是：Sum（坐标*强度值）/ Sum（坐标）
将每个群集移动到新的质心。

这种方法几乎可以肯定地提供更好的结果，并且您可以获得每个群集的质量，这可能有助于识别脚趾。

（同样，您已经预先指定了群集的数量。使用群集，您必须以这样或那样的方式指定密度：在这种情况下选择适合的群集数量，或者选择群集半径并查看结束的群集数量后者的一个例子是平均移位。）

对于缺乏实施细节或其他细节感到抱歉。我会对此进行编码，但我有一个截止日期。如果下周没有其他任何工作让我知道，我会试一试。

#9楼

看来你可以使用jetxee的算法作弊。他发现前三个脚趾很好，你应该能够猜到第四个脚趾的位置。

#10楼

我不确定这回答了这个问题，但看起来你可以找到没有邻居的n个最高峰。

这是要点。请注意，它是在Ruby中，但这个想法应该是清楚的。

require 'pp'NUM_PEAKS = 5
NEIGHBOR_DISTANCE = 1data = [[1,2,3,4,5],[2,6,4,4,6],[3,6,7,4,3],]def tuples(matrix)tuples = []matrix.each_with_index { |row, ri|row.each_with_index { |value, ci|tuples << [value, ri, ci]}}tuples
enddef neighbor?(t1, t2, distance = 1)[1,2].each { |axis|return false if (t1[axis] - t2[axis]).abs > distance}true
end# convert the matrix into a sorted list of tuples (value, row, col), highest peaks first
sorted = tuples(data).sort_by { |tuple| tuple.first }.reverse# the list of peaks that don't have neighbors
non_neighboring_peaks = []sorted.each { |candidate|# always take the highest peakif non_neighboring_peaks.empty?non_neighboring_peaks << candidateputs "took the first peak: #{candidate}"else# check that this candidate doesn't have any accepted neighborsis_ok = truenon_neighboring_peaks.each { |accepted|if neighbor?(candidate, accepted, NEIGHBOR_DISTANCE)is_ok = falsebreakend}if is_oknon_neighboring_peaks << candidateputs "took #{candidate}"elseputs "denied #{candidate}"endend
}pp non_neighboring_peaks

#11楼

这是我在为大型望远镜做类似事情时使用的另一种方法：

1）搜索最高像素。一旦你有了这个，搜索它最适合2x2（可能最大化2x2总和），或在4x4的子区域内以最高像素为中心进行2d高斯拟合。

然后将您发现的2x2像素设置为峰值中心周围的零（或3x3）

回到1）并重复直到最高峰值低于噪音阈值，或者你有所需的所有脚趾

#12楼

如果你能够创建一些训练数据，那么尝试使用神经网络可能是值得的......但是这需要手工注释许多样本。

#13楼

粗略的轮廓......

您可能希望使用连通组件算法来隔离每个爪区域。维基对此有一个很好的描述（有一些代码）： http ： //en.wikipedia.org/wiki/Connected_Component_Labeling

你必须决定是否使用4或8连通性。就个人而言，对于大多数问题，我更喜欢6连通性。无论如何，一旦你将每个“爪印”分离为一个连接区域，它应该很容易迭代整个区域并找到最大值。一旦找到了最大值，就可以迭代地放大该区域，直到达到预定的阈值，以便将其识别为给定的“脚趾”。

这里的一个微妙问题是，一旦你开始使用计算机视觉技术来识别某些东西作为右/左/前/后爪并且你开始查看单个脚趾，你必须开始考虑旋转，倾斜和翻译。这是通过分析所谓的“时刻”来实现的。在视觉应用中需要考虑几个不同的时刻：

中心矩：平移不变的归一化矩：缩放和平移不变的时刻：平移，缩放和旋转不变

有关时刻的更多信息可以通过在wiki上搜索“图像时刻”找到。

#14楼

也许你可以使用像Gaussian Mixture Models这样的东西。这是一个用于执行GMM的Python程序包（只是进行了Google搜索） http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

#15楼

物理学家的解决方案：
定义由其位置X_i标识的5个爪标记，并用随机位置初始化它们。定义一些能量函数，将一些奖励用于爪子位置的标记位置，并对标记的重叠进行一些惩罚; 让我们说：

E(X_i;S)=-Sum_i(S(X_i))+alfa*Sum_ij (|X_i-Xj|<=2*sqrt(2)?1:0)

（ S(X_i)是围绕X_i 2x2平方的平均力， alfa是通过实验达到峰值的参数）

现在是时候做一些Metropolis-Hastings魔术了：
1.选择随机标记并沿随机方向将其移动一个像素。
2.计算dE，此移动引起的能量差异。
3.从0-1获得一个统一的随机数并将其称为r。
4.如果dE<0或exp(-beta*dE)>r ，接受移动并转到1; 如果没有，撤消移动并转到1。
这应该重复，直到标记会聚到爪子。 Beta控制扫描以优化权衡，因此它也应该通过实验进行优化; 它也可以随着模拟时间（模拟退火）不断增加。

#16楼

有趣的问题。我想尝试的解决方案如下。

应用低通滤波器，例如使用2D高斯掩模进行卷积。这将为您提供一堆（可能，但不一定是浮点）值。
使用每个爪垫（或脚趾）的已知近似半径执行2D非最大抑制。

这应该给你最大的位置，而没有多个候选人在一起。只是为了澄清，步骤1中掩模的半径也应该类似于步骤2中使用的半径。这个半径可以是可选择的，或者兽医可以事先明确地测量它（它将随着年龄/品种/等而变化）。

建议的一些解决方案（平均移位，神经网络等）可能会在某种程度上起作用，但过于复杂，可能并不理想。

#17楼

物理学家已经在一定程度上研究了这个问题。 ROOT有一个很好的实现。查看TSpectrum类（特别是针对您的情况的TSpectrum2 ）及其文档。

参考文献：

M.Morhac等人：用于多维重合伽马射线谱的背景消除方法。核仪器和物理研究方法A 401（1997）113-132。
M.Morhac等：高效的一维和二维金去卷积及其在伽马射线光谱分解中的应用。核仪器和物理研究方法A 401（1997）385-408。
M.Morhac等人：在多维重合伽马射线光谱中识别峰。研究物理中的核仪器和方法A 443（2000），108-125。

......对于那些无权订阅NIM的人：

Spectrum.doc
SpectrumDec.ps.gz
SpectrumSrc.ps.gz
SpectrumBck.ps.gz

#18楼

感谢原始数据。我在火车上，这是我已经到达的（我的站点即将到来）。我用regexps按摩了你的txt文件，然后把它放到一个带有一些javascript的html页面中进行可视化。我在这里分享它是因为有些人，比如我自己，可能会发现它比python更容易被黑客攻击。

我认为一个好的方法是比例和旋转不变，我的下一步将是调查高斯的混合物。（每个爪垫是高斯的中心）。

    <html>
<head><script type="text/javascript" src="http://vis.stanford.edu/protovis/protovis-r3.2.js"></script> <script type="text/javascript">var heatmap = [[[0,0,0,0,0,0,0,4,4,0,0,0,0],
[0,0,0,0,0,7,14,22,18,7,0,0,0],
[0,0,0,0,11,40,65,43,18,7,0,0,0],
[0,0,0,0,14,61,72,32,7,4,11,14,4],
[0,7,14,11,7,22,25,11,4,14,65,72,14],
[4,29,79,54,14,7,4,11,18,29,79,83,18],
[0,18,54,32,18,43,36,29,61,76,25,18,4],
[0,4,7,7,25,90,79,36,79,90,22,0,0],
[0,0,0,0,11,47,40,14,29,36,7,0,0],
[0,0,0,0,4,7,7,4,4,4,0,0,0]
],[
[0,0,0,4,4,0,0,0,0,0,0,0,0],
[0,0,11,18,18,7,0,0,0,0,0,0,0],
[0,4,29,47,29,7,0,4,4,0,0,0,0],
[0,0,11,29,29,7,7,22,25,7,0,0,0],
[0,0,0,4,4,4,14,61,83,22,0,0,0],
[4,7,4,4,4,4,14,32,25,7,0,0,0],
[4,11,7,14,25,25,47,79,32,4,0,0,0],
[0,4,4,22,58,40,29,86,36,4,0,0,0],
[0,0,0,7,18,14,7,18,7,0,0,0,0],
[0,0,0,0,4,4,0,0,0,0,0,0,0],
],[
[0,0,0,4,11,11,7,4,0,0,0,0,0],
[0,0,0,4,22,36,32,22,11,4,0,0,0],
[4,11,7,4,11,29,54,50,22,4,0,0,0],
[11,58,43,11,4,11,25,22,11,11,18,7,0],
[11,50,43,18,11,4,4,7,18,61,86,29,4],
[0,11,18,54,58,25,32,50,32,47,54,14,0],
[0,0,14,72,76,40,86,101,32,11,7,4,0],
[0,0,4,22,22,18,47,65,18,0,0,0,0],
[0,0,0,0,4,4,7,11,4,0,0,0,0],
],[
[0,0,0,0,4,4,4,0,0,0,0,0,0],
[0,0,0,4,14,14,18,7,0,0,0,0,0],
[0,0,0,4,14,40,54,22,4,0,0,0,0],
[0,7,11,4,11,32,36,11,0,0,0,0,0],
[4,29,36,11,4,7,7,4,4,0,0,0,0],
[4,25,32,18,7,4,4,4,14,7,0,0,0],
[0,7,36,58,29,14,22,14,18,11,0,0,0],
[0,11,50,68,32,40,61,18,4,4,0,0,0],
[0,4,11,18,18,43,32,7,0,0,0,0,0],
[0,0,0,0,4,7,4,0,0,0,0,0,0],
],[
[0,0,0,0,0,0,4,7,4,0,0,0,0],
[0,0,0,0,4,18,25,32,25,7,0,0,0],
[0,0,0,4,18,65,68,29,11,0,0,0,0],
[0,4,4,4,18,65,54,18,4,7,14,11,0],
[4,22,36,14,4,14,11,7,7,29,79,47,7],
[7,54,76,36,18,14,11,36,40,32,72,36,4],
[4,11,18,18,61,79,36,54,97,40,14,7,0],
[0,0,0,11,58,101,40,47,108,50,7,0,0],
[0,0,0,4,11,25,7,11,22,11,0,0,0],
[0,0,0,0,0,4,0,0,0,0,0,0,0],
],[
[0,0,4,7,4,0,0,0,0,0,0,0,0],
[0,0,11,22,14,4,0,4,0,0,0,0,0],
[0,0,7,18,14,4,4,14,18,4,0,0,0],
[0,4,0,4,4,0,4,32,54,18,0,0,0],
[4,11,7,4,7,7,18,29,22,4,0,0,0],
[7,18,7,22,40,25,50,76,25,4,0,0,0],
[0,4,4,22,61,32,25,54,18,0,0,0,0],
[0,0,0,4,11,7,4,11,4,0,0,0,0],
],[
[0,0,0,0,7,14,11,4,0,0,0,0,0],
[0,0,0,4,18,43,50,32,14,4,0,0,0],
[0,4,11,4,7,29,61,65,43,11,0,0,0],
[4,18,54,25,7,11,32,40,25,7,11,4,0],
[4,36,86,40,11,7,7,7,7,25,58,25,4],
[0,7,18,25,65,40,18,25,22,22,47,18,0],
[0,0,4,32,79,47,43,86,54,11,7,4,0],
[0,0,0,14,32,14,25,61,40,7,0,0,0],
[0,0,0,0,4,4,4,11,7,0,0,0,0],
],[
[0,0,0,0,4,7,11,4,0,0,0,0,0],
[0,4,4,0,4,11,18,11,0,0,0,0,0],
[4,11,11,4,0,4,4,4,0,0,0,0,0],
[4,18,14,7,4,0,0,4,7,7,0,0,0],
[0,7,18,29,14,11,11,7,18,18,4,0,0],
[0,11,43,50,29,43,40,11,4,4,0,0,0],
[0,4,18,25,22,54,40,7,0,0,0,0,0],
[0,0,4,4,4,11,7,0,0,0,0,0,0],
],[
[0,0,0,0,0,7,7,7,7,0,0,0,0],
[0,0,0,0,7,32,32,18,4,0,0,0,0],
[0,0,0,0,11,54,40,14,4,4,22,11,0],
[0,7,14,11,4,14,11,4,4,25,94,50,7],
[4,25,65,43,11,7,4,7,22,25,54,36,7],
[0,7,25,22,29,58,32,25,72,61,14,7,0],
[0,0,4,4,40,115,68,29,83,72,11,0,0],
[0,0,0,0,11,29,18,7,18,14,4,0,0],
[0,0,0,0,0,4,0,0,0,0,0,0,0],
]
];
</script>
</head>
<body><script type="text/javascript+protovis">    for (var a=0; a < heatmap.length; a++) {var w = heatmap[a][0].length,h = heatmap[a].length;
var vis = new pv.Panel().width(w * 6).height(h * 6).strokeStyle("#aaa").lineWidth(4).antialias(true);
vis.add(pv.Image).imageWidth(w).imageHeight(h).image(pv.Scale.linear().domain(0, 99, 100).range("#000", "#fff", '#ff0a0a').by(function(i, j) heatmap[a][j][i]));
vis.render();
}
</script></body>
</html>

#19楼

我使用局部最大滤波器检测到峰值。以下是您的第一个4爪数据集的结果：

我还在9个爪子的第二个数据集上运行它，它也运行良好。

这是你如何做到的：

import numpy as np
from scipy.ndimage.filters import maximum_filter
from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
import matplotlib.pyplot as pp#for some reason I had to reshape. Numpy ignored the shape header.
paws_data = np.loadtxt("paws.txt").reshape(4,11,14)#getting a list of images
paws = [p.squeeze() for p in np.vsplit(paws_data,4)]def detect_peaks(image):"""Takes an image and detect the peaks usingthe local maximum filter.Returns a boolean mask of the peaks (i.e. 1 whenthe pixel's value is the neighborhood maximum, 0 otherwise)"""# define an 8-connected neighborhoodneighborhood = generate_binary_structure(2,2)#apply the local maximum filter; all pixel of maximal value #in their neighborhood are set to 1local_max = maximum_filter(image, footprint=neighborhood)==image#local_max is a mask that contains the peaks we are #looking for, but also the background.#In order to isolate the peaks we must remove the background from the mask.#we create the mask of the backgroundbackground = (image==0)#a little technicality: we must erode the background in order to #successfully subtract it form local_max, otherwise a line will #appear along the background border (artifact of the local maximum filter)eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)#we obtain the final mask, containing only peaks, #by removing the background from the local_max mask (xor operation)detected_peaks = local_max ^ eroded_backgroundreturn detected_peaks#applying the detection and plotting results
for i, paw in enumerate(paws):detected_peaks = detect_peaks(paw)pp.subplot(4,2,(2*i+1))pp.imshow(paw)pp.subplot(4,2,(2*i+2) )pp.imshow(detected_peaks)pp.show()

您需要做的就是在掩码上使用scipy.ndimage.measurements.label来标记所有不同的对象。然后你就可以单独玩它们了。

请注意 ，该方法效果很好，因为背景不嘈杂。如果是，你会在后台检测到一堆其他不需要的峰值。另一个重要因素是邻里的规模。如果峰值大小发生变化，则需要对其进行调整（应保持大致成比例）。

#20楼

使用持久同源性来分析您的数据集我得到以下结果（点击放大）：

这是本SO答案中描述的峰值检测方法的2D版本。上图仅显示了按持久性排序的0维持久同源类。

我使用scipy.misc.imresize（）将原始数据集升级了2倍。但请注意，我确实将四个爪子视为一个数据集; 将它分成四个可以使问题更容易。

方法。 这背后的想法非常简单：考虑将每个像素分配给它的级别的函数的函数图。它看起来像这样：

现在考虑高度为255的水位不断下降到较低水位。在当地最大岛屿上弹出（出生）。在马鞍点两个岛屿合并; 我们认为下岛要合并到更高的岛屿（死亡）。所谓的持久性图（第0维同源类，我们的岛屿）描绘了所有岛屿的死亡率 - 超过出生值：

那么岛屿的持续存在就是出生和死亡水平之间的差异; 点与灰色主对角线的垂直距离。该图通过减少持久性来标记岛屿。

第一张照片显示了这些岛屿的出生地点。该方法不仅给出局部最大值，而且通过上述持久性量化它们的“重要性”。然后，人们会过滤掉持久性太低的所有岛屿。但是，在您的示例中，每个岛屿（即每个局部最大值）都是您要寻找的峰值。

Python代码可以在这里找到。

#21楼

只是想告诉你们有一个很好的选择，用python找到图像中的局部最大值。

from skimage.feature import peak_local_max

或者对于skimage 0.8.0

from skimage.feature.peak import peak_local_max

http://scikit-image.org/docs/0.8.0/api/skimage.feature.peak.html

#22楼

天文学和宇宙学界提供了多种广泛的软件 - 这在历史和现在都是一个重要的研究领域。

如果你不是天文学家，请不要惊慌 - 有些人很容易在野外使用。例如，您可以使用astropy / photutils：

https://photutils.readthedocs.io/en/stable/detection.html#local-peak-detection

[在这里重复他们的简短示例代码似乎有点粗鲁。]

下面给出了可能感兴趣的不完整且略有偏见的技术/包/链接列表 - 请在评论中添加更多内容，我将根据需要更新此答案。当然，存在精确度与计算资源的权衡。 [老实说，在这样的单一答案中提供代码示例太多了，所以我不确定这个答案是否会成功。]

源提取器https://www.astromatic.net/software/sextractor

MultiNest https://github.com/farhanferoz/MultiNest [+ pyMultiNest]

ASKAP / EMU源寻找挑战： https ：//arxiv.org/abs/1509.03931

您还可以搜索普朗克和/或WMAP源提取挑战。

...