二维阵列中的峰值检测

 

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

我制作了每个爪子的二维数组，每个爪子随着时间的推移装载了每个传感器的最大值。 这里有一个爪子的例子，我用Excel绘制了我想检测的区域。 这些传感器周围有2 x 2个盒子，带有当地最大值，总和最大。

替代文字

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

替代文字

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

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

此外，为了方便起见，我还采用了2x2，任何更高级的解决方案都是值得欢迎的，但我只是一名人类运动科学家，所以我既不是真正的程序员，也不是数学家，所以请保持简单。

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

结果

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

更具体地说，它不能识别第四趾的小峰。 这显然是固有的，因为循环从最低值开始往下看，而没有考虑到这一点。

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

替代文字

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

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

替代文字

更新：

我已经为任何感兴趣的人建立了一个博客，并且我已经用所有原始尺寸设置了SkyDrive。 所以对于任何需要更多数据的人来说：给你更多的权力

新的更新：

所以在帮助我了解了有关爪子检测和爪子分类的问题之后，我终于能够检查每只爪子的脚趾检测了！ 事实证明，除了我自己的例子中的那种尺寸的爪子，它的效果并不好。 事后看来，如此任意选择2x2是我自己的错。

这里有一个很好的例子，它出现了问题：一个钉子被认为是一个脚趾，而'脚跟'很宽，它会被识别两次！

替代文字

爪子太大，所以采取2×2大小没有重叠，导致一些脚趾被检测两次。 反过来说，小狗通常不会找到第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平方的总和。 我认为你需要所有这些，所以可能会有一些重叠。 我使用切片从原始二维数组中切出第一/最后一列和最后一行，然后将它们重叠在一起并计算总和。

为了更好地理解它，对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]])

现在想象你将它们叠加在一起，并将元素加在相同的位置。 这些总和将与在同一位置左上角的2×2方格中的总和完全相同： 

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

当你有超过2×2的正方形时，你可以使用max来查找最大值，或者sort ，或者sorted来查找峰值。

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

笔记

我允许2×2的方块重叠。 编辑过的版本会过滤掉其中的一些，这样只有非重叠的正方形出现在结果中。

选择手指（一个想法）

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

我注意到，如果前手指几乎保持一个完美的圆形，后手指应该在该圆内。 另外，前手指或多或少地等间距。 我们可能会尝试使用这些启发式属性来检测手指。

伪代码： 

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均值聚类。 它没有那么坏。

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

    • 然后迭代直到收敛：

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

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

    （同样，你已经指定了簇的数量，在聚类中你必须以这样或那样的方式指定密度：在这种情况下选择簇的数量，或者选择一个簇半径，看看你结束了多少后者的一个例子是mean-shift。）

    对于缺乏实现细节或其他细节感到抱歉。 我会编码，但我有一个截止日期。 如果下周没有其他工作让我知道，我会给它一个镜头。

    链接地址: http://www.djcxy.com/p/4533.html

    上一篇: Peak detection in a 2D array

    下一篇: How does Python's super() work with multiple inheritance?