二维阵列中的峰值检测

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

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

alt text

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

alt text

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

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

另外，为了方便起见，我选择了2x2，欢迎任何更先进的解决方案，但我只是一个人类运动科学家，所以我既不是真正的程序员也不是数学家，所以请保持“简单”。

这是一个可以加载#0的版本

搜索结果

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

更具体地说，它无法识别第四个脚趾的小峰。这显然是固有的事实，即循环自上而下地看向最低值，而不考虑它在哪里。

有没有人知道如何调整@jextee的算法，以便它也能找到第四个脚趾？

alt text

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

这张照片显示了它们是如何在空间上分布在盘子上的。

alt text

更新时间：

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

新增更新：

因此，在我得到关于爪检测和爪分选的问题的帮助后，我终于能够检查每个爪子的脚趾检测！事实证明，除了像我自己的例子中那样大小的爪子之外，它在任何地方都不起作用。事后看来，如此武断地选择2x2是我自己的错。

这里有一个很好的例子，说明它出了什么问题：指甲被识别为脚趾，而“脚后跟”太宽了，它被识别了两次！

alt text

爪子太大了，所以采取2x2大小没有重叠，导致一些脚趾被检测到两次。另一方面，在小型犬中，它经常找不到第五个脚趾，我怀疑这是由于2x2区域太大造成的。

在在我所有的测量中尝试当前的解决方案之后，我得出了一个惊人的结论，几乎所有我的小狗都没有找到第五个脚趾，而在超过50%的大型狗的影响中，它会找到更多！

所以很明显我需要改变它。我自己的猜测是将neighborhood的大小更改为小型狗的小尺寸和大型狗的大尺寸。但是generate_binary_structure不让我改变数组的大小。

因此，我希望其他人有更好的建议来定位脚趾，也许有脚趾面积与爪子大小的比例？

134823 次浏览
小开

只是我脑海中的几个想法：

  • 取扫描的梯度（导数），看看是否可以消除错误调用
  • 取局部最大值的最大值

你可能还想看看OpenCV，它有一个相当不错的Python API，可能有一些你会觉得有用的函数。

小开

也许一个简单的方法就足够了：在你的平面上建立一个所有2x2方块的列表，按它们的总和排序（降序）。

首先，在你的“爪子列表”中选择价值最高的方块。然后，迭代选择4个与之前找到的任何方块都不相交的次优方块。

小开

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

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

小开

如果你一步一步地进行：你首先找到全局最大值，如果需要处理周围的点，给出它们的值，然后将找到的区域设置为零，然后重复下一个。

小开

解决方案

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

from scipy import *from operator import itemgetter
n = 5  # how many fingers are we looking for
d = 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 squares
def 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_overlaps
pairs2 = drop_overlapping(sorted(pairs, key=itemgetter(1), reverse=True))
# Take the first n with the heighest values
positions = pairs2[:n]
# Print results
print 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) ; aarray([[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:]; sumsarray([[ 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例。

小开

好吧，这里有一些简单而不是非常有效的代码，但对于这种大小的数据集来说，它是好的。

import numpy as npgrid = 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个左右的初始设置。单独运行每个方向，并丢弃两个或多个脚趾最终位于同一像素的任何结果。我会再考虑一下，但是这种事情仍在图像处理中研究-没有正确的答案！

稍微复杂一点的想法：（加权）K均值聚类。没那么糟糕。

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

然后迭代直到收敛：

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

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

（同样，你已经在前面指定了集群的数量。使用集群，你必须以某种方式指定密度：要么选择集群的数量，在这种情况下合适，要么选择集群半径，看看你最终得到多少。后者的一个例子是均值漂移。）

很抱歉缺少实现细节或其他细节。我会把它编码起来，但我有一个截止日期。如果下周没有其他工作，请告诉我，我会试一试。

小开

看起来你可以用jetxee的算法作弊，他发现前三个脚趾很好，你应该能猜出第四个脚趾在哪里。

小开

我不确定这是否回答了这个问题，但似乎你可以只寻找没有邻居的n个最高峰。

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

require 'pp'
NUM_PEAKS = 5NEIGHBOR_DISTANCE = 1
data = [[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]}}tuplesend
def neighbor?(t1, t2, distance = 1)[1,2].each { |axis|return false if (t1[axis] - t2[axis]).abs > distance}trueend
# convert the matrix into a sorted list of tuples (value, row, col), highest peaks firstsorted = tuples(data).sort_by { |tuple| tuple.first }.reverse
# the list of peaks that don't have neighborsnon_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
小开

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

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

然后将峰值中心周围的2x2像素设置为零（或者可能是3x3）

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

小开

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

小开

一个粗略的轮廓…

您可能希望使用连通分支算法来隔离每个爪子区域。wiki在这里对此有一个不错的描述（带有一些代码）：http://en.wikipedia.org/wiki/Connected_Component_Labeling

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

这里的一个微妙问题是，一旦你开始使用计算机视觉技术来识别某物为右/左/前/后爪子，并开始观察单个脚趾，你就必须开始考虑旋转、倾斜和平移。这是通过分析所谓的“时刻”来完成的。在视觉应用中有几个不同的时刻需要考虑：

中心矩：平移不变量归一化矩：缩放和平移不变量hu矩：平移、缩放和旋转不变量

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

小开

也许你可以使用高斯混合模型之类的东西。这是一个用于做GMM的Python包（刚刚做了谷歌搜索）http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

小开

物理学家的解决方案：
定义5个由其位置X_i标识的爪子标记，并使用随机位置初始化它们。定义一些能量函数，将标记在爪子位置的一些奖励与标记重叠的一些惩罚结合起来；让我们说：

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是实验达到峰值的参数）

现在是时候做一些大都会黑斯廷斯魔术了：
1.选择随机标记并向随机方向移动一个像素。
2.计算dE，这一步造成的能量差。
3.从0-1中获取一个统一的随机数，并将其称为r。
4.如果dE<0或exp(-beta*dE)>r，接受移动并转到1；如果不是，撤消移动并转到1。
这应该重复，直到标记将收敛到爪子。Beta控制扫描以优化权衡，所以它也应该在实验上优化；它也可以随着模拟（模拟退火）的时间不断增加。

小开

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

  1. 应用低通滤波器，例如使用2D高斯掩模进行卷积。这将为您提供一堆（可能，但不一定是浮点）值。

  2. 使用每个爪垫（或脚趾）的已知近似半径执行2D非最大抑制。

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

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

小开

物理学家已经对这个问题进行了深入的研究。在ROOT中有一个很好的实现。看看TSpectrum类（尤其是TSpectrum2对于你的情况）和它们的留档。

参考文献：

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

…对于那些无法订阅NIM的人：

小开

谢谢你的原始数据。我在火车上，这是我所得到的（我的站点即将到来）。我用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>

alt文本

小开
最佳答案

我使用局部最大滤波检测到峰值。这是你的第一个4只爪子数据集的结果：峰值检测结果

我还在第二个数据集上运行了9个爪子和它也起作用了。

以下是你如何做到这一点：

import numpy as npfrom scipy.ndimage.filters import maximum_filterfrom scipy.ndimage.morphology import generate_binary_structure, binary_erosionimport 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 imagespaws = [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_background
return detected_peaks

#applying the detection and plotting resultsfor 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来标记所有不同的对象。然后你就可以单独玩它们了。

说明表示该方法效果良好，因为背景没有噪音。如果是，你会在背景中检测到一堆其他不需要的峰值。另一个重要因素是街区的大小。如果峰值大小发生变化，你需要调整它（应该保持大致成比例）。

小开

我相信你现在有足够的时间继续下去，但是我忍不住建议使用k均值聚类方法。k均值是一种无监督聚类算法，它会把你的数据（在任何维度上——我碰巧是在3D上做的）并将其排列成k个具有不同边界的集群。这里很好，因为你确切地知道这些犬科动物（应该）有多少脚趾。

此外，它是在非常好的Slipy中实现的（http://docs.scipy.org/doc/scipy/reference/cluster.vq.html）。

以下是它可以在空间上解析3D集群的示例：在此处输入图像描述

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

小开

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

结果

这是本所以回答中描述的峰值检测方法的2D版本。上图简单地显示了按持久性排序的0维持久同源类。

我确实使用scipy.misc.imresize（）将原始数据集提升了2倍。然而，请注意，我确实将四个爪子视为一个数据集；将其拆分为四个会使问题更容易。

<强>方法论。这背后的想法非常简单：考虑为每个像素分配其级别的函数图。它看起来像这样：

3D函数图

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

持久化图

岛屿的持久性是出生级和死亡级之间的差异；一个点到灰色主对角线的垂直距离。该图通过减少持久性来标记岛屿。

第一张图片显示了岛屿的诞生位置。这种方法不仅给出了局部最大值，还通过上述持久性量化了它们的“显著性”。然后，人们会过滤掉所有持久性太低的岛屿。然而，在你的例子中，每个岛屿（即每个局部最大值）都是你寻找的峰值。

可以找到Python代码这里。

小开

只是想告诉你们，有一个不错的选择可以使用python在图像中找到本地maxima：

from skimage.feature import peak_local_max

对于skimage0.8.0：

from skimage.feature.peak import peak_local_max

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


提交答案