python - skimage安装 - scikit-image安装




二维阵列中的峰值检测 (14)

我正在帮助兽医诊所测量狗爪下的压力。 我使用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不会让我改变数组的大小。

因此,我希望其他人对脚趾定位有更好的建议,也许脚趾面积与脚掌尺寸一致?


数据文件: 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]])

现在想象你将它们叠加在一起,并将元素加在相同的位置。 这些总和将与在同一位置左上角的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个情况。


I am not sure this answers the question, but it seems like you can just look for the n highest peaks that don't have neighbors.

Here is the gist. Note that it's in Ruby, but the idea should be clear.

require 'pp'

NUM_PEAKS = 5
NEIGHBOR_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]
    }
  }
  tuples
end

def neighbor?(t1, t2, distance = 1)
  [1,2].each { |axis|
    return false if (t1[axis] - t2[axis]).abs > distance
  }
  true
end

# convert the matrix into a sorted list of tuples (value, row, col), highest peaks first
sorted = tuples(data).sort_by { |tuple| tuple.first }.reverse

# the list of peaks that don't have neighbors
non_neighboring_peaks = []

sorted.each { |candidate|
  # always take the highest peak
  if non_neighboring_peaks.empty?
    non_neighboring_peaks << candidate
    puts "took the first peak: #{candidate}"
  else
    # check that this candidate doesn't have any accepted neighbors
    is_ok = true
    non_neighboring_peaks.each { |accepted|
      if neighbor?(candidate, accepted, NEIGHBOR_DISTANCE)
        is_ok = false
        break
      end
    }
    if is_ok
      non_neighboring_peaks << candidate
      puts "took #{candidate}"
    else
      puts "denied #{candidate}"
    end
  end
}

pp non_neighboring_peaks

It seems you can cheat a bit using jetxee's algorithm. He is finding the first three toes fine, and you should be able to guess where the fourth is based off that.


Maybe a naive approach is sufficient here: Build a list of all 2x2 squares on your plane, order them by their sum (in descending order).

First, select the highest-valued square into your "paw list". Then, iteratively pick 4 of the next-best squares that don't intersect with any of the previously found squares.


一个粗略的轮廓...

你可能想使用连接组件算法来隔离每个爪子区域。 wiki有一个体面的描述(与一些代码)在这里: http://en.wikipedia.org/wiki/Connected_Component_Labeling : http://en.wikipedia.org/wiki/Connected_Component_Labeling

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

一个微妙的问题是,只要你开始使用计算机视觉技术来识别右/左/前/后爪的东西,并且你开始看单个脚趾,就必须考虑旋转,倾斜和平移。 这是通过分析所谓的“时刻”来完成的。 视觉应用中需要考虑几个不同的时刻:

中心矩:平移不变标准化矩:缩放和平移不变胡矩:平移,缩放和旋转不变

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


下面介绍另一种在大型望远镜上做类似事情时使用的方法:

1)搜索最高像素。 一旦你有了,找出最适合2x2(也许最大化2x2和)的方法,或者在以最高像素为中心的4x4子区域内进行2d高斯拟合。

然后在峰值中心周围设置你发现的2×2像素为零(或可能是3×3)

回到1)并重复,直到最高峰落在噪声阈值以下,或者您拥有所有需要的脚趾


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

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

我使用scipy.misc.imresize()将原始数据集放大了2倍。 但是,请注意,我确实将四个爪子视为一个数据集; 将它分成四个会使问题更容易。

方法。 这背后的想法很简单:考虑函数的函数图表,为每个像素指定其级别。 它看起来像这样:

现在考虑一个高度为255的水位,该水位不断下降到较低水平。 在当地最大的岛屿出现(出生)。 在鞍点两个岛屿合并; 我们认为较低的岛屿将被合并到较高的岛屿(死亡)。 所谓的持久性图(属于第0维的同源类,我们的岛屿)描述了所有岛屿的死亡超过的出生值:

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

第一张照片显示了这些岛屿的出生地点。 这种方法不仅给出了局部极大值,而且还通过上述持久性量化了它们的“显着性”。 然后,人们会过滤掉所有持续时间太短的岛屿。 然而,在你的例子中,每个岛屿(即每个本地最大值)都是你寻找的高峰。



如果您能够创建一些训练数据,那么尝试使用神经网络可能是值得的...但是这需要手动注释许多样本。


感谢原始数据。 我在火车上,这是我所得到的(我的停下来了)。 我用正则表达式处理了你的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>


我相信你现在已经够了,但我不禁建议使用k-means聚类方法。 k-means是一种无监督聚类算法,可以将数据(任意维数 - 我恰巧以3D方式进行处理)获取数据,并将其排列为具有不同边界的k个群集。 这里很好,因为你知道这些犬(应该)有多少脚趾。

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

以下是可以在空间上解析3D群集的一个示例:

你想要做的是有点不同(2D和包括压力值),但我仍然认为你可以给它一个镜头。


物理学家已经深入研究了这个问题。 ROOT有一个很好的实现。 看看TSpectrum类(特别是针对您的案例的TSpectrum2 )以及它们的文档。

参考文献:

  1. M.Morhac等:多维重合伽马射线谱的背景消除方法。 Nuclear Instruments and Methods in Physics Research A 401(1997)113-132。
  2. M.Morhac等:高效的一维和二维金解卷积及其在伽玛射线光谱分解中的应用。 Nuclear Instruments and Methods in Physics Research A 401(1997)385-408。
  3. M.Morhac等:在多维重合伽马射线谱中的峰的鉴定。 Nuclear Instruments and Methods in Research Physics A 443(2000),108-125。

...对于那些无法获得NIM订阅的用户:


这是一个图像注册问题 。 总的策略是:

  • 有一个已知的例子,或某种事先的数据。
  • 将您的数据放在示例中,或者将示例适用于您的数据。
  • 如果您的数据首先大致对齐,它会有所帮助。

这是一个粗略和准备的方法 ,“可能有用的最愚蠢的东西”:

  • 大致在您期望的地方从五个脚趾坐标开始。
  • 随着每一个,迭代爬到山顶。 即给定当前位置,如果其值大于当前像素,则移动到最大相邻像素。 当你的脚趾坐标停止移动时停止。

为了抵消定位问题,您可以对基本方向(北,东北等)进行8次左右的初始设置。 分别运行每一个,并丢弃任何结果,其中两个或更多个脚趾以相同像素结束。 我会更多地考虑这一点,但这种事情仍在图像处理中进行研究 - 没有正确的答案!

稍微更复杂的想法:(加权)K均值聚类。 它没有那么坏。

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

然后迭代直到收敛:

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

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

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

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


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

这里有另外一个想法:如果你知道高压点的典型大小,你可以先用相同大小的高斯函数对它进行卷积来平滑你的图像。 这可能会让您处理更简单的图像。







image-processing