在一个二维数组中进行峰值检测。

1001
我正在帮助一家兽医诊所测量狗爪下的压力。我使用Python进行数据分析,但现在我遇到了一个问题,即如何将爪子划分为(解剖学上的)子区域。
我创建了一个二维数组,其中包含每个爪子的最大值,这些最大值是由爪子随时间加载的传感器产生的。以下是一个示例,我使用Excel绘制了我想要“检测”的区域。这些区域是围绕具有局部极大值的传感器的2x2方框,它们的总和最大。

alt text

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

alt text

这其中哪一个最大值是我想要的,告诉Python最好的方式是什么?
注意:2x2的方块不能重叠,因为它们必须是分开的。
我把2x2当作一种便利方法,欢迎任何更高级的解决方案,但我只是一个人类运动科学家,既不是真正的程序员也不是数学家,所以请保持“简单”。
这是一个可以使用`np.loadtxt`加载的版本

结果

所以我尝试了 @jextee 的解决方案(请参见下面的结果)。正如您所看到的,它在前脚上运作得非常好,但对于后腿来说效果不太好。

更具体地说,它无法识别作为第四趾的小峰。这显然是由于该循环从上向下看最低值,并未考虑其位置造成的固有问题。

是否有人知道如何调整 @jextee 的算法,使其能够找到第四趾?

alt text

因为我还没有处理其他试验,所以无法提供其他样本。但之前给出的数据是每只爪子的平均值。这个文件是一个数组,按照它们与板子接触的顺序列出了9只爪子的最大数据。

这张图显示了它们在板子上如何空间分布。

alt text

更新:

我为所有感兴趣的人建立了一个博客,并且我已经设置了一个包含所有原始测量数据的OneDrive。因此,对于任何请求更多数据的人:祝你好运!


新的更新:

在我关于 爪子检测爪子分类 的问题得到帮助后,我终于能够检测每只爪子的脚趾!结果发现,除了与我自己示例中的爪子尺寸相似的爪子外,它在其他情况下的效果并不好。事后想想,选择2x2的尺寸是我的任性。

以下是一个错误的例子:一根指甲被识别为脚趾,而“脚跟”太宽,被识别了两次!

alt text

爪子太大了,所以采用2x2的尺寸且不重叠,会导致有些脚趾被检测到两次。反过来,在小狗身上通常无法找到第五个脚趾,我怀疑是因为2x2区域太大的原因。

对所有我的测量值尝试当前解决方案后,令人惊讶的结论是几乎所有的小狗都没有找到第五个脚趾,而在大狗的超过50%的情况下会找到更多!

所以显然我需要进行改变。我自己的猜测是将neighborhood的大小对于小狗设定得更小,对于大狗设定得更大。但是generate_binary_structure不允许我更改数组的大小。

有人对于定位脚趾有更好的建议吗?也许能够使脚趾区域与爪子大小成比例变化?


1
我理解逗号是小数点而不是值分隔符,对吗? - MattH
1
是的,它们是逗号。@Christian,我试图将其放入一个易于阅读的文件中,但即使这样也失败了 :( - Ivo Flipse
5
由于我正在进行可行性研究,因此任何方法都可以。所以我正在寻找尽可能多的定义压力的方法,包括亚区域。此外,我需要能够区分“大脚趾”和“小脚趾”两侧,以估计方向。但由于这之前没有做过,所以我们无法预测会发现什么 :-) - Ivo Flipse
3
@Ron:这项研究的目标之一是看这个系统适用于什么体型/重量的狗,所以是的,虽然这只狗大约有20公斤。我还有一些明显较小(和更大)的狗,并预计我不会能够对真正小的狗做同样的事情。 - Ivo Flipse
3
@frank,爪子的尺寸是随时间测量的,因此有第三个维度。然而,它们不会移动位置(相对而言),所以我更关注爪子在二维空间中的位置。在此基础上,第三个维度是额外得到的。 - Ivo Flipse
显示剩余5条评论
21个回答

390

我使用了本地最大值过滤器来检测峰值。以下是在第一个包含4只爪子的数据集上的结果: Peaks detection result

我还在包含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 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """

    # define an 8-connected neighborhood
    neighborhood = generate_binary_structure(2,2)

    #apply the local maximum filter; all pixel of maximal value 
    #in their neighborhood are set to 1
    local_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 background
    background = (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 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函数进行标记以便单独处理。

请注意,该方法能够奏效的原因是背景不会出现噪点。如果有噪点的话,你将会检测到许多其他不需要的波峰。另一个重要因素是邻域的大小。如果波峰的大小发生变化,你需要相应地调整它(应该保持大致成比例)。


3
比起(eroded_background ^ local_peaks),有一个更简单的解决方案。只需执行(foreground & local peaks)即可。 - Ryan Soklaski

60

解决方案

数据文件: 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) ; 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查找最大值,或者使用sortsorted查找峰值。
为了记住峰值的位置,我将每个值(总和)与其在平铺数组中的序数位置配对(参见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, radius
        check if the selected finger is inside of the circle
        check 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种组合,每个爪子需要评估3960种情况,因为每个爪子有5种选择后指的方式。


鉴于您的结果列表,他将不得不手动过滤掉这些无用选项...选择前四个最高的结果将为他提供构建一个最大值为6.8的2x2方格的四个可能性。 - Johannes Charra
1
我试过了,对于前爪似乎可行,但对于后爪效果不佳。看来我们必须尝试一些知道该去哪里寻找的东西。 - Ivo Flipse
我看到问题所在了。我会思考如何识别最佳的“星座”峰值来选择。你觉得采用“四连一旁”或“环形四连一内”的方法怎么样? - sastanin
正如我第二张图片所示(这里有所有爪子的链接),如果您检查每行和每列的最大值,则所有峰值都会被标记,因此也许我们不仅可以按从上到下排序的方式遍历列表,而且还可以检查哪些最大值是最高的,同时没有邻居(忽略靠近最大值的一切)。甚至可以查看每行和每列的2x2总和哪个最大。 - Ivo Flipse
如果我们使用一些启发式方法来确定最有可能的两个最高脚趾的候选者,也许可以根据形状来确定后脚趾,这样就可以减少组合的数量。此外,从其他建议中观察高斯滤波器的使用情况,或许可以增加您的建议的有效性。 - Ivo Flipse

40

这是一个图像配准问题。一般的策略是:

  • 有一个已知的示例,或者对数据有某种先验知识。
  • 将数据拟合到示例上,或将示例拟合到数据上。
  • 如果你的数据在第一次大致对齐,会有帮助。

以下是一种粗略但可行的方法:

  • 从五个脚趾坐标开始,大致放置在你期望的位置上。
  • 对于每个脚趾坐标,迭代地向山顶攀登。例如,给定当前位置,移动到相邻的像素中最大的一个,如果其值大于当前像素,则继续移动。直到脚趾坐标不再移动为止。

为了解决方向性问题,你可以有约8个初始设置,分别代表基本方向(北、东北等)。单独运行每个设置,并且舍弃任何两个或以上脚趾在同一像素处的结果。我会进一步考虑这个问题,但这种方法在图像处理中仍在研究中-没有正确答案!

稍微更复杂的想法:(加权)K均值聚类。并不那么困难。

  • 从五个脚趾坐标开始,但现在这些是“聚类中心”。

然后迭代直到收敛:

  • 将每个像素分配给最近的聚类(为每个聚类创建一个列表)。
  • 计算每个聚类的重心。对于每个聚类,这是:Sum(坐标*强度值)/Sum(坐标)
  • 将每个聚类移动到新的重心。

这种方法几乎肯定会提供更好的结果,并且你可以得到每个聚类的质量,这可能有助于识别脚趾。

(再次强调,您提前指定了聚类数量。在聚类中,您必须以某种方式指定密度:要么选择适合本案例的聚类数量,要么选择聚类半径并查看您最终得到了多少聚类。后者的一个例子是均值漂移。)

很抱歉没有实现细节或其他具体信息。我会编写代码,但我有截止日期。如果下周之前没有任何进展,请告诉我,我会尝试一下。


32

使用持续同调来分析您的数据集,我得到了以下结果(点击放大):

Result

这是描述在SO answer中的峰值检测方法的二维版本。上图仅显示按持续性排序的0维持久同调类。
我使用scipy.misc.imresize()将原始数据集放大了2倍。然而,请注意,我将四只脚视为一个数据集;将其分割成四个将使问题更容易。
方法论。 这背后的思想相当简单:考虑将每个像素赋予其级别的函数的函数图。它看起来像这样:

3D function graph

现在考虑一条高度为255的水平线,不断下降到更低的水平。在局部最大值处会出现小岛(诞生)。在鞍点处,两个小岛会合并;我们认为较低的小岛与较高的小岛合并(死亡)。所谓的持久图(0维同调类的小岛)描述了所有小岛的死亡-诞生值:

Persistence diagram

岛屿的持久性是其出生水平和死亡水平之间的差异;即点到灰色主对角线的垂直距离。图中通过递减的持久性标记了岛屿。

第一张图片展示了岛屿的出生位置。该方法不仅提供局部极大值,还通过上述持久性量化它们的“重要性”。然后可以过滤掉所有持久性过低的岛屿。但是,在您的示例中,每个岛屿(即每个局部最大值)都是您寻找的峰值。

Python代码可以在此处找到。


1
我在C++中实现了相同的算法,比答案中链接的Python实现快约45倍。 C++实现可在此处找到:https://gitlab.com/balping/find-peaks - balping

21
这个问题已经被物理学家深入研究过。在 ROOT 中有一个很好的实现。看看 TSpectrum 类(特别是 TSpectrum2 适用于你的情况)和它们的文档。
参考文献:
1. M.Morhac 等人:多维符合伽马射线谱背景消除方法。《核仪器和物理研究方法 A》401 (1997) 113-132。 2. M.Morhac 等人:高效的一维和二维 Gold 去卷积及其在 γ 射线谱分解中的应用。《核仪器和物理研究方法 A》401 (1997) 385-408。 3. M.Morhac 等人:多维符合伽马射线谱峰识别。《核仪器和物理研究方法 A》443(2000), 108-125。
... 对于那些没有 NIM 订阅的人:

19
我不禁建议使用k-means聚类方法。k-means是一种无监督的聚类算法,它可以将您的数据(在任意维度上 - 我碰巧在3D中进行)分成k个具有明确边界的簇。这里很好的原因是您知道这些犬科动物(应该)有多少只脚趾。
此外,它在Scipy中实现得非常好(http://docs.scipy.org/doc/scipy/reference/cluster.vq.html)。
以下是一个示例,展示了它如何在空间上解决3D簇: enter image description here 虽然您想要做的事情有点不同(2D并包括压力值),但我仍然认为您可以尝试一下。

16

这里有一个想法:计算图像的(离散)拉普拉斯算子。我预计它在极值点处会更加突出,并且更加剧烈,相对于原始图像。因此,可能更容易找到极值点。

这里还有另一个想法:如果你知道高压点的典型大小,可以先通过与同样大小的高斯核卷积平滑图像。这样可能会给您更简单的图像来处理。


15

以下是我脑海里想到的一些点子:

  • 获取扫描图像的梯度(导数),看能否消除虚假信号
  • 获取局部最大值的最大值

您可能还想看一下 OpenCV,它有一个相当不错的Python API,并且可能有一些您会发现有用的函数。


使用梯度,您的意思是我应该计算斜坡的陡峭程度,一旦超过某个值,我就知道有“山峰”了?我尝试过这样做,但有些脚趾只有非常低的山峰(1.2 N/cm),与其他一些脚趾相比(8 N/cm),如何处理具有非常低梯度的山峰? - Ivo Flipse
2
过去对我有用的方法是,如果我不能直接使用梯度,就要查看梯度和最大值,例如,如果梯度是局部极值,并且我在局部最大值处,则我就到达了一个感兴趣的点。 - ChrisC

13
我用正则表达式处理了你的文本文件,并将其放入一个带有一些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 text


12

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接