二维阵列中的峰检测

我正在帮助兽医诊所测量狗爪下的压力。我使用 Python 进行数据分析,现在我被困在试图将爪子分成(解剖)子区域。

我制作了每个爪子的 2D 数组,其中包含爪子随时间推移已加载的每个传感器的最大值。这是一个爪子的示例,我使用 Excel 绘制了要 “检测” 的区域。这些是传感器周围具有最大最大值的 2 x 2 框,它们的总和最大。

替代文字

因此,我尝试了一些实验,并决定只寻找每一列和每一行的最大值(由于爪子的形状而无法朝一个方向看)。这似乎可以很好地 “检测” 到各个脚趾的位置,但是它也标记了相邻的传感器。

替代文字

那么,告诉 Python 我想要的最大最大值是什么?

注意:2x2 的正方形不能重叠,因为它们必须是单独的脚趾!

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

是可以通过np.loadtxt加载版本


结果

因此,我尝试了 @jextee 的解决方案(请参见下面的结果)。如您所见,它在前爪上非常有效,但在后腿上效果较差。

更具体地说,它无法识别出第四脚趾的小峰。显然,这是回路从上向下朝着最低值而不考虑这在哪里的事实所固有的。

谁会知道如何调整 @jextee 的算法,以便它也能够找到第四个脚趾?

替代文字

由于我尚未处理其他任何试验,因此无法提供任何其他样品。但是我之前提供的数据是每只爪子的平均值。该文件是一个数组,其中最大 9 爪的数据按它们与板接触的顺序排列。

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

替代文字

更新:

我已经为有兴趣的任何人建立了博客,为 SkyDrive 设置了所有原始测量值。因此,对于任何需要更多数据的人:给您更多的权力!


新更新:

因此,在获得有关爪子检测爪子分类的问题的帮助后,我终于能够检查每个爪子的脚趾检测!事实证明,除了爪子大小像我自己的示例中的爪子一样,它在任何情况下都无法正常工作。事后看来,如此随意地选择 2x2 是我自己的错。

这是一个出问题的好例子:指甲被识别为脚趾,而 “脚跟” 是如此之宽,被识别两次!

替代文字

脚掌太大,因此采用 2x2 大小且没有重叠的脚掌会使两次脚趾被检测到两次。相反,在小型犬中,它通常找不到第 5 个脚趾,我怀疑这是 2x2 区域太大造成的。

尝试了所有解决方案的最新解决方案后,我得出了一个惊人的结论:几乎对我所有的小型犬来说,它都找不到第五个脚趾,而在大型犬的 50%以上的撞击中,它会发现更多!

所以很明显我需要更改它。我自己的猜测是将neighborhood的大小更改为小型犬较小,大型犬较大。但是generate_binary_structure不允许我更改数组的大小。

因此,我希望其他人对脚趾的定位有更好的建议,也许脚趾的面积与爪子的大小成正比?

答案

我使用局部最大滤波器检测到峰值。这是第一个 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 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标记所有不同的对象。这样您就可以分别与他们一起玩了。

请注意 ,该方法效果很好,因为背景不嘈杂。如果是这样,您将在背景中检测到许多其他不需要的峰。另一个重要因素是邻里的大小。如果峰大小发生变化,则需要对其进行调整(应保持大致成比例)。

数据文件: 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 种组合,乘以 5 种选择后指的方式,因此每只爪子需要评估 3960 种情况。

这是图像配准问题 。总体策略是:

  • 有一个已知的例子,或某种先验的数据。
  • 使数据适合示例,或使示例适合数据。
  • 如果您的数据首先是大致对齐的,则将有帮助。

这是一个粗略而现成的方法 ,“可能会起作用的最愚蠢的事情”:

  • 从五个脚趾坐标开始,大致在您期望的位置。
  • 与每个人一起,反复攀登到山顶。即给定当前位置,如果其值大于当前像素,则移动到最大相邻像素。当脚趾坐标停止移动时停止。

要解决方向问题,您可以为基本方向(北,东北等)设置 8 个左右的初始设置。单独运行每个,并丢弃两个或多个脚趾最终位于同一像素的任何结果。我会再考虑一下,但是在图像处理中仍在研究这种东西 - 没有正确的答案!

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

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

然后迭代直到收敛:

  • 将每个像素分配给最近的聚类(只需为每个聚类列出一个列表)。
  • 计算每个群集的质心。对于每个群集,这是:总和(坐标 * 强度值)/ 总和(坐标)
  • 将每个群集移动到新的质心。

这种方法几乎可以肯定会带来更好的结果,而且您会得到每个簇的质量,这可能有助于识别脚趾。

(同样,您已经预先指定了群集数。使用群集时,您必须以一种或另一种方式指定密度:选择适合这种情况的群集数,或者选择一个群集半径,然后查看要终止的群集数后者的一个例子是均值漂移 。)

抱歉,缺少实施细节或其他细节。我会对此进行编码,但是我有最后期限。如果下周没有其他工作,请告诉我,我会尝试一下。

物理学家已经对该问题进行了深入研究。在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 的用户:

使用持久同源性分析您的数据集,我得到以下结果(点击放大):

结果

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

我使用 scipy.misc.imresize()将原始数据集按比例放大了 2 倍。但是,请注意,我确实将四个爪子视为一个数据集。将其分为四个部分将使问题更容易解决。

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

3D功能图

现在考虑一个高度为 255 的水位,该水位连续下降到较低的水位。在局部最大的岛上弹出(出生)。在马鞍点,两个岛屿合并;我们认为较低的岛屿将合并到较高的岛屿(死亡)。所谓的持久性图(属于零维度同质性类,我们的岛屿)描述了所有岛屿的死亡人数与出生人数的关系:

持续图

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

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

Python 代码可以在这里找到。

这是一个想法:计算图像的(离散)拉普拉斯算子。我希望它在最大值上(负)很大,比原始图像更具戏剧性。因此,最大值可能更容易找到。

这是另一个想法:如果您知道高压点的典型大小,则可以先通过使用相同大小的高斯将其卷积来平滑图像。这可以使您处理更简单的图像。

我的脑海中只有几个想法:

  • 进行扫描的梯度(导数),看看是否消除了错误的调用
  • 取局部最大值的最大值

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

我确定您现在有足够的工作要做,但是我不禁建议使用 k-means 聚类方法。 k-means 是一种无监督的聚类算法,它将获取您的数据(任意数量的维度 - 我碰巧是在 3D 中进行此操作)并将其排列到具有不同边界的 k 个聚类中。这里很好,因为您确切知道这些犬(应该)有多少个脚趾。

此外,它是在 Scipy 中实现的,这真的很好( http://docs.scipy.org/doc/scipy/reference/cluster.vq.html )。

这是如何在空间上解析 3D 群集的示例: 在此处输入图片说明

您想要做的有点不同(2D 并包含压力值),但是我仍然认为您可以尝试一下。

感谢您的原始数据。我在火车上,这是我所能到达的最大距离(我的停靠站即将到来)。我用正则表达式对您的 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>

替代文字

物理学家的解决方案:
定义 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)是在 2×2 的周围方均值力X_ialfa是通过实验达到顶峰的参数)

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