二维阵列中的峰值检测

Question

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

我为每个paw制作了一个2D数组，它包含了随着时间推移paw加载的每个传感器的最大值。这里有一个单爪的例子，我用Excel来绘制我想“检测”的区域。这些是传感器周围的2×2盒子，带有局部最大值，它们的总和最大。

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

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

注意：2x2方块不能重叠，因为它们必须是分开的脚趾！

另外，我把2x2作为一个方便，任何更先进的解决方案都是受欢迎的，但我只是一个人类运动科学家，所以我既不是一个真正的程序员，也不是一个数学家，所以请保持它'简单'。

这是一张version that can be loaded with ^{}

---

结果

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

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

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

因为我还没有做过其他的试验，所以我不能提供其他的样品。但我之前给出的数据是每个爪子的平均值。这个文件是一个数组，最大数据为9个爪子，按它们与板接触的顺序排列。

这张图片显示了它们是如何在空间上分布在板块上的。

更新：

I have set up a blog for anyone interested和I have setup a SkyDrive with all the raw measurements.所以对任何要求更多数据的人来说：给你更多的力量！

---

新更新：

所以在我得到关于paw detection和paw sorting的问题的帮助后，我终于可以检查每个爪子的脚趾检测了！结果，它在任何东西上都没有那么好用，只是爪子大小和我自己例子中的一样。事后看来，这是我自己的错，选择2x2如此武断。

这里有一个很好的例子说明它哪里出错了：指甲被识别为脚趾，而‘脚跟’太宽了，它被识别两次！

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

在trying the current solution on all my measurements之后，我得出了一个惊人的结论：对于我所有的小狗来说，它都没有找到第五个脚趾，而在超过50%的撞击中，对于大狗来说，它会找到更多！

很明显我需要改变它。我自己的猜测是把neighborhood的大小改成小狗的小，大狗的大。但是generate_binary_structure不允许我更改数组的大小。

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

Answer 1

这是一个image registration problem。总体战略是：

• 有一个已知的例子，或者数据上的某种先验。
• 根据示例调整数据，或根据数据调整示例。
• 如果您的数据最初是大致对齐的，那么它会有帮助。

这里有一个粗略且准备就绪的方法，“可能最愚蠢的事情”：

• 从五个脚趾坐标开始，大致在你期望的位置。
• 每一个都要反复爬到山顶。i、 e.给定当前位置，如果其值大于当前像素，则移动到最大相邻像素。当你的脚趾坐标停止移动时停止。

为了解决方向问题，可以为基本方向（北、东北等）设置8个左右的初始设置。分别运行每一个并丢弃两个或多个脚趾在同一像素处结束的任何结果。我会再想一想，但这类事情还在图像处理领域研究——没有正确的答案！

稍微复杂一点的想法：（加权）K-均值聚类。还不错。

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

然后迭代直到收敛：

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

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

（同样，您已经指定了前面的集群数量。对于集群，您必须以一种或另一种方式指定密度：要么选择集群的数量（在本例中是适当的），要么选择集群的半径并查看最终的数量。后者的一个例子是mean-shift）

很抱歉，缺少实现细节或其他细节。我会把它编码，但我有最后期限。如果到下个星期为止没有别的办法，请告诉我，我会试一试的。

Answer 2

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

我还在第二个数据集上运行了它，它包含9个paws和it worked as well。

以下是您的操作方法：

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来标记所有不同的对象。然后你就可以单独和他们玩了。

注意该方法工作良好，因为背景没有噪音。如果是的话，你会在背景中发现一堆不需要的峰。另一个重要因素是邻域的大小。如果峰值大小发生变化，则需要对其进行调整（应保持大致成比例）。

Answer 3

溶液

数据文件：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"

Output没有重叠的正方形。似乎与您的示例中选择的区域相同。

一些评论

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

要更好地理解它，请对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, 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例。