二维数组中的峰值检测

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

我做了每个爪子的二维数组,每个爪子随着时间的推移装载的每个传感器的最大值。 下面是一个爪子的例子,我用Excel来绘制我想要检测的区域。 这些都是与本地最大值的传感器周围2乘2箱,总共有最大的总和。

替代文字

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

替代文字

那么告诉Python哪个最好的是我想要的呢?

注意:2×2方块不能重叠,因为它们必须是分开的脚趾!

另外我还带了2×2作为一个方便,更高级的解决scheme是值得欢迎的,但我只是一个人类运动科学家,所以我既不是一个真正的程序员或math家,所以请保持简单。

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


结果

所以我尝试了@ jextee的解决scheme(见下面的结果)。 正如你所看到的,它在前爪上非常有效,但后腿的效果不太好。

更具体地说,它不能识别第四脚趾的小峰。 这显然是固有的事实,即循环看起来低于最低值,而没有考虑到这一点。

有谁会知道如何调整@ jextee的algorithm,以便它能够find第四个脚趾吗?

替代文字

由于我还没有处理任何其他的审判,我不能提供任何其他样本。 但是我之前给出的数据是每个爪子的平均值。 这个文件是一个数组,其最大数据量为9个爪子,它们与盘子接触的顺序是这样的。

这幅图像显示了它们是如何在板上空间分布的。

替代文字

更新:

我已经为所有感兴趣的人build立了博客,并且我已经用所有的原始尺寸设置了SkyDrive。 所以要求更多数据的人:给你更多的权力!


更新:

所以在我的帮助后,我得到了关于爪子检测和爪子sorting的问题 ,终于能够检查每个爪子的脚趾检测! 事实certificate,除了我自己的例子中的那种尺寸的爪子之外,它的效果并不好。 事后看来,如此任意select2x2是我自己的错。

这里有一个很好的例子:一个钉子被认为是一个脚趾,而“脚跟”是如此之宽,它会被认出两次!

替代文字

爪子太大,所以采取2×2大小没有重叠,导致一些脚趾被检测到两次。 反过来说,小狗往往没有find第五个脚趾,我怀疑是由于2×2面积太大造成的。

在对我所有的测量结果进行了现在的解决之后,我得出了惊人的结论,几乎所有的小型犬都没有find第五只脚趾,而在超过50%的大型狗的影响中,它会发现更多!

很显然,我需要改变它。 我自己的猜测是把neighborhood的大小改变成小狗的大小,大狗的大小。 但是, generate_binary_structure不会让我改变数组的大小。

因此,我希望其他人对脚趾定位有更好的build议,也许脚趾的尺寸与脚掌尺寸有关?

我使用一个局部最大值filter来检测峰值。 这是你的第一个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 (ie 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" 

输出不重叠的正方形。 看来你的例子中select了相同的区域。

一些意见

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

要更好地理解它,对3x3arrays进行成像:

 >>> 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]]) 

现在想象一下,将它们叠加在一起,然后在相同的位置上求和元素。 这些总和在与左上angular相同的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来find峰值。

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

笔记

我允许2×2的方块重叠。 编辑版本将其中的一部分过滤掉,结果只有非重叠的正方形出现。

select手指(一个想法)

另一个问题是如何从所有的峰值中select可能的手指。 我有一个想法,可能会或可能不会工作。 我没有时间来执行它,所以只是伪代码。

我注意到,如果前手指几乎保持完美的一圈,那么后手指应该在该圈内。 而且,正面的手指或多或less是等距的。 我们可以尝试使用这些启发式属性来检测手指。

伪代码:

 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种select后手指的方式,因此对每个爪子评估3960个情况。

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

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

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

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

为了解决方向问题,你可以有8个左右的基本方向(北,东北等)的初始设置。 单独运行每一个,并丢弃任何结果,其中两个或更多个脚趾以相同的像素结束。 我会多想一些,但是这种东西还在image processing中被研究 – 没有正确的答案!

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

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

然后迭代,直到收敛:

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

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

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

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

这个问题已经被物理学家深入研究了。 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的用户:

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

只是我的头顶上的几个想法:

  • 采取扫描的梯度(派生),看看是否消除了错误的电话
  • 采取当地最大值的最大值

您可能还想看看OpenCV ,它有一个相当不错的Python API,并可能有一些function,你会发现有用的。

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

这里有另一个想法:如果你知道高压点的典型大小,你可以先用一个相同大小的高斯来使图像平滑。 这可能会给你更简单的image processing。

物理学家的解决scheme:
定义由它们的位置X_i标识的5个爪标记,并用随机位置初始化它们。 定义一些能量函数,结合某些标记在爪子位置上的位置奖励以及对标记重叠的一些惩罚; 让我们说:

 E(X_i;S)=-Sum_i(S(X_i))+alfa*Sum_ij (|X_i-Xj|<=2*sqrt(2)?1:0) 

S(X_i)S(X_i)周围2×2平方的平均力, alfa是实验达到峰值的参数)

现在时间做一些大都会 – 黑斯廷斯的魔法:
1.select随机标记并随机移动一个像素。
2.计算dE,这个移动引起的能量差异。
3.从0-1得到一个统一的随机数,并称之为r。
4.如果dE<0exp(-beta*dE)>r ,则接受移动并转到1; 如果没有,撤消移动并转到1。
这应该重复,直到标记将收敛到爪子。 Beta控制扫描以优化折衷,因此它也应该被实验优化; 随着模拟时间的增加(模拟退火)也可以不断增加。

感谢您的原始数据。 我在火车上,而且就我所知(我停下来了)。 我使用正则expression式来处理你的txt文件,并用一些javascript将其放到一个html页面中以便显示。 我在这里分享它,因为像我自己一样,可能会发现它比python更容易hackable。

我认为一个好的方法将是规模和旋转不变,我的下一步将是调查高斯混合物。 (每个掌垫都是高斯的中心)。

  <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> 

替代文字

另一种我在做大型望远镜类似工作时使用的方法:

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

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

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

如果能够创build一些训练数据,那么尝试使用neural network可能是值得的…但是这需要手动注释许多样本。

粗略的轮廓

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

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

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

中心矩:平移不变归一化矩:尺度和平移不变胡矩:平移,尺度和旋转不variables

更多关于时刻的信息可以通过在wiki上search“图像时刻”来find。

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

 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] 

我基本上只是用左上angular的位置和每个2×2平方和的数组sorting并按总和sorting。 然后我把2×2的方块与最高的总和竞争,把它放在best数组中,并删除所有其他使用这个刚刚移除的2×2方块的任何部分的2×2方块。

它似乎工作正常,除了最后一个爪子(在第一张图片的最右边的总和最小的那个),事实certificate还有另外两个合格的2×2正方形具有较大的总和(并且它们具有等于彼此)。 其中之一仍然是从你的2×2平方中select一个方格,但是另一个方向是在左边。 幸运的是,幸运的是,我们看到的是更多地select你想要的,但是这可能需要一些其他的想法来得到你实际上想要的东西。

看来你可以使用jetxee的algorithm作弊。 他发现前三个脚趾都没问题,你应该能猜出第四个脚趾的位置。

也许你可以使用像高斯混合模型的东西。 下面是一个用来做GMM的Python软件包(刚做了一个Googlesearch) http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

有趣的问题。 我会尝试的解决scheme如下。

  1. 应用一个低通滤波器,例如用二维高斯遮罩进行卷积。 这会给你一堆(可能,但不一定是浮点)值。

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

这应该给你最大的职位,没有多个候选人紧密在一起。 为了澄清,步骤1中的面罩的半径也应该类似于步骤2中使用的半径。该半径可以是可select的,或者兽医可以事先明确地测量它(它将随着年龄/品种/等而变化)。

一些解决scheme(平均移位,neural network等)可能会在一定程度上起作用,但是过于复杂并且可能不理想。

我相信你现在已经够了,但我不禁build议使用k-means聚类方法。 k-means是一种无监督的聚类algorithm,它将把你的数据(在任何维数 – 我碰巧做到这一点在三维),并安排到不同的边界k个集群。 这里很好,因为你知道这些犬(应该)有多less脚趾。

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

这是一个例子,它可以做什么来空间parsing3D群集: 在这里输入图像描述

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

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

结果

这是在这个答案中描述的峰值检测方法的2D版本。 上图简单展示了通过持久性sorting的0维持久同源类。

I did upscale the original dataset by a factor of 2 using scipy.misc.imresize(). However, note that I did consider the four paws as one dataset; splitting it into four would make the problem easier.

Methodology. The idea behind this quite simple: Consider the function graph of the function that assigns each pixel its level. 它看起来像这样:

3D function graph

Now consider a water level at height 255 that continuously descents to lower levels. At local maxima islands pop up (birth). At saddle points two islands merge; we consider the lower island to be merged to the higher island (death). The so-called persistence diagram (of the 0-th dimensional homology classes, our islands) depicts death- over birth-values of all islands:

Persistence diagram

The persistence of an island is then the difference between the birth- and death-level; the vertical distance of a dot to the grey main diagonal. The figure labels the islands by decreasing persistence.

The very first picture shows the locations of births of the islands. This method not only gives the local maxima but also quantifies their "significance" by the above mentioned persistence. One would then filter out all islands with a too low persistence. However, in your example every island (ie, every local maximum) is a peak you look for.

Maybe a naive approach is sufficient here: Build a list of all 2×2 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.

What if you proceed step by step: you first locate the global maximum, process if needed the surrounding points given their value, then set the found region to zero, and repeat for the next one.

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