二维阵列中的峰值检测

我用激光探测到了这些峰 局部最大滤波器 . 以下是您的第一个4爪数据集的结果:

我还对第二个数据集9爪和 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

注意 邻里 . 如果峰值大小发生变化,则需要对其进行调整(应保持大致成比例)。

解决方案

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

一些评论

为了更好地理解它,请对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个病例。

这是一个 image registration problem

• 有一个已知的例子,或者某种 先前的 在数据上。
• 将数据与示例相匹配,或将示例与数据相匹配。
• 如果你的数据是 首先对齐。

这里有一个粗略的方法

• 从五个脚趾的坐标开始,大致在你期望的地方。

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

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

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

然后迭代直至收敛:

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

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

mean-shift

很抱歉,缺少实现细节或其他细节。我本想把它编出来的,但我还有最后期限。如果下星期没有别的办法奏效,请告诉我,我会试试看的。

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

我确实使用scipy.misc.imresize调整大小(). 但是,请注意,我确实将四个paws视为一个数据集;将其拆分为四个将使问题更容易解决。

方法论。

现在考虑一个高度为255的水位,它不断下降到较低的水位。在当地的马克西马群岛弹出(出生)。在鞍点两个岛屿合并;我们认为较低的岛屿合并到较高的岛屿(死亡)。所谓的持久性图(第0维同源类,我们的岛屿)描述了所有岛屿的死亡-出生值:

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

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

Python代码可以找到 here .

ROOT . 看看这个 TSpectrum 课程(尤其是 TSpectrum2

参考文献:

2. M.Morhac等人:有效的一维和二维黄金反卷积及其在伽马射线谱分解中的应用。物理学研究中的核仪器和方法a401(1997)385-408。
3. M.Morhac等人:多维符合伽马射线谱峰的识别。研究物理学中的核仪器和方法A 443(2000),108-125。

…对于那些无法订阅NIM的用户:

• Spectrum.doc
• SpectrumDec.ps.gz
• SpectrumSrc.ps.gz
• SpectrumBck.ps.gz

这里有一个想法:计算图像的(离散)拉普拉斯函数。我希望它是(消极的)最大,在一种方式,这是更戏剧性的比原来的形象大。因此,maxima可能更容易找到。

我脑子里有几个想法:

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

你也许还想看看 OpenCV ,它有一个相当不错的pythonapi,可能有一些你会发现有用的函数。

另外,它是在Scipy中实现的,这非常好( http://docs.scipy.org/doc/scipy/reference/cluster.vq.html ).

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

谢谢你的原始数据。我在火车上,这是我到的最远的地方(我的站快到了)。我用regexps处理了你的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>

物理学家的解决方案:
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) Xïi alfa 是一个要在实验中达到峰值的参数)

现在是时候做一些大都会黑斯廷斯魔术了:

2计算dE,这一移动引起的能量差。

dE<0 或 exp(-beta*dE)>r
这应该重复,直到标记将收敛到爪子。Beta控制扫描以优化折衷,因此它也应该通过实验进行优化;它也可以随着模拟时间(模拟退火)而不断增加。

我只想告诉你们有个很好的选择 maxima 在使用python的图像中:

from skimage.feature import peak_local_max

还是为了撇渣 0.8.0

from skimage.feature.peak import peak_local_max

http://scikit-image.org/docs/0.8.0/api/skimage.feature.peak.html

下面是另一种方法,我在对大型望远镜进行类似操作时使用:

1) 搜索最高像素。 一旦你有了这个,搜索2x2的最佳拟合(也许最大化2x2和),或者在4x4的子区域内做一个2d高斯拟合,以最高像素为中心。

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

大概的轮廓。。。

您可能希望使用连接组件算法来隔离每个paw区域。wiki对此有很好的描述(带有一些代码): http://en.wikipedia.org/wiki/Connected_Component_Labeling

你必须决定是使用4还是8连接性。就我个人而言,对于大多数问题,我更喜欢6连接性。不管怎样,一旦你把每个“爪印”分离出来作为一个连通的区域,遍历这个区域并找到最大值就足够容易了。一旦找到了最大值,就可以迭代地放大区域,直到达到预定的阈值,以便将其识别为给定的“脚趾”。

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

归一化矩:缩放和平移不变性 hu矩:平移、缩放和旋转不变

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

也许你可以使用高斯混合模型。下面是一个用于执行GMMs的Python包(刚刚进行了Google搜索) http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

有趣的问题。我想尝试的解决方案如下。

1. 应用低通滤波器,例如用二维高斯掩模卷积。这将为您提供一组(可能,但不一定是浮点)值。

2. 使用每个爪垫(或脚趾)的已知近似半径执行2D非最大抑制。

这应该给你最大的位置,而不是有多个候选人是紧密在一起。只是澄清一下,步骤1中遮罩的半径也应该与步骤2中使用的半径相似。这个半径可以选择,或者兽医可以事先明确测量它(它会随着年龄/品种等而变化)。

建议的一些解决方案(均值漂移、神经网络等)可能会在某种程度上起作用,但过于复杂,可能并不理想。

好吧,这里有一些简单而不是非常有效的代码,但是对于这样大小的数据集来说,这是很好的。

import numpy as np
grid = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0],
              [0,0,0,0,0,0,0,0,0.4,0.4,0.4,0,0,0],
              [0,0,0,0,0.4,1.4,1.4,1.8,0.7,0,0,0,0,0],
              [0,0,0,0,0.4,1.4,4,5.4,2.2,0.4,0,0,0,0],
              [0,0,0.7,1.1,0.4,1.1,3.2,3.6,1.1,0,0,0,0,0],
              [0,0.4,2.9,3.6,1.1,0.4,0.7,0.7,0.4,0.4,0,0,0,0],
              [0,0.4,2.5,3.2,1.8,0.7,0.4,0.4,0.4,1.4,0.7,0,0,0],
              [0,0,0.7,3.6,5.8,2.9,1.4,2.2,1.4,1.8,1.1,0,0,0],
              [0,0,1.1,5,6.8,3.2,4,6.1,1.8,0.4,0.4,0,0,0],
              [0,0,0.4,1.1,1.8,1.8,4.3,3.2,0.7,0,0,0,0,0],
              [0,0,0,0,0,0.4,0.7,0.4,0,0,0,0,0,0]])

arr = []
for i in xrange(grid.shape[0] - 1):
    for j in xrange(grid.shape[1] - 1):
        tot = grid[i][j] + grid[i+1][j] + grid[i][j+1] + grid[i+1][j+1]
        arr.append([(i,j),tot])

best = []

arr.sort(key = lambda x: x[1])

for i in xrange(5):
    best.append(arr.pop())
    badpos = set([(best[-1][0][0]+x,best[-1][0][1]+y)
                  for x in [-1,0,1] for y in [-1,0,1] if x != 0 or y != 0])
    for j in xrange(len(arr)-1,-1,-1):
        if arr[j][0] in badpos:
            arr.pop(j)


for item in best:
    print grid[item[0][0]:item[0][0]+2,item[0][1]:item[0][1]+2]

我基本上只是用左上角的位置和每个2x2平方的和做一个数组,然后按和排序。然后我把2x2平方和的最大值从争用中取出来,放到 best 数组,并删除所有其他2x2正方形,使用任何部分,这只是删除2x2正方形。

除了最后一只爪子(第一张图片中最右边和最小的一只爪子)之外,它似乎工作得很好,结果发现还有另外两个合格的2x2平方和更大(它们的和彼此相等)。其中一个仍然是从2x2正方形中选择一个正方形,但另一个在左侧关闭。幸运的是,幸运的是,我们可以从你想要的东西中选择更多,但这可能需要一些其他的想法来得到你一直想要的东西。

我不确定这是否回答了这个问题,但似乎你可以只寻找没有邻居的n个最高峰。

Here is the gist. 注意,它是用Ruby编写的,但是这个想法应该很清楚。

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

https://photutils.readthedocs.io/en/stable/detection.html#local-peak-detection

[在这里重复他们的简短示例代码似乎有点粗鲁。]

下面给出了一个不完整的和有点偏颇的可能感兴趣的技术/包/链接列表-请在评论中添加更多内容,我会根据需要更新这个答案。当然,精度和计算资源之间有一个权衡。[老实说,像这样的一个答案中有太多的代码示例,所以我不确定这个答案会不会飞起来。]

源提取器 https://www.astromatic.net/software/sextractor

多巢 https://github.com/farhanferoz/MultiNest [+pyMultiNest]

ASKAP/EMU寻源挑战: https://arxiv.org/abs/1509.03931

您还可以搜索普朗克和/或WMAP源提取挑战。

如果您一步一步地进行:首先定位全局最大值,如果需要,处理给定值的周围点,然后将找到的区域设置为零,然后重复下一个区域。