numpy:大量线段/点的快速有规律的平均间隔 [英] numpy: fast regularly-spaced average for large numbers of line segments / points

问题描述

我沿着一维线有许多(〜一百万个)不规则间隔的点P.这些标记线段，使得如果点是{0，x_a，x_b，x_c，x_d，...}，则线段从0-> x_a，x_a-> x_b，x_b-> x_c，x_c- > x_d等.对于每个片段，我都有一个y值，我希望将其解释为颜色的深度.我需要将这条线绘制为图像，但是可能只有(比如说)1000个像素可以代表该线的整个长度.当然，这些像素对应于沿着线的规则间隔，例如在0..X1，X1..X2，X2..X3等处，其中X1，X2，X3是规则间隔的.要计算每个像素的颜色，我需要取落在规则间隔的像素边界内的所有y值的平均值，并按该间隔内段的长度加权.可能还有一些像素在P中不包含任何值，它们只是采用了跨整个像素的线段所定义的颜色值.

I have many (~ 1 million) irregularly spaced points, P, along a 1D line. These mark segments of the line, such that if the points are {0, x_a, x_b, x_c, x_d, ...}, the segments go from 0->x_a, x_a->x_b, x_b->x_c, x_c->x_d, etc. I also have a y-value for each segment, which I wish to interpret as a depth of colour. I need to plot this line as an image, but there might only be (say) 1000 pixels available to represent the entire length of the line. These pixels, of course, correspond to regularly spaced intervals along the line, say at 0..X1, X1..X2, X2..X3, etc, where X1, X2, X3 are regularly spaced. To work out the colour for each pixel, I need to take an average of all the y-values that fall within the regularly-spaced pixel boundaries, weighted by the length of the segment that falls in that interval. There are also likely to be pixels which do not contain any value in P, and which simply take the colour value defined by the segment that passes across the entire pixel.

似乎在图像分析中可能需要做很多事情.那么，该操作有一个名称吗?在numpy中，什么是最快的方法来计算这样的一组规则的y平均值间隔?我想这有点像插值法，只是我不想只取两个周围点的平均值，而是取规则间隔内所有点的加权平均值(加上一点重叠).

This seems like something that probably needs to be done a lot in image analysis. So is there a name for this operation, and what's the fastest way in numpy to calculate such a regularly-spaced set of average y-values? It's a bit like interpolation, I guess, only I don't want to take the mean of just the two surrounding points, but a weighted average of all points within a regular interval (plus a bit of overlap).

假设一条水平线上有5个段，它们由[0，1.1，2.2，2.3，2.8，4]分隔(即，该行从0到4).假设每个线段采用任意阴影值，例如，我们可以使用5个阴影值[0,0.88,0.55,0.11,0.44]-其中0为黑色，而1为白色.然后，如果我想使用4个像素进行绘制，则需要从0 ... 1、1 ... 2等创建4个值，并希望计算结果为每个值返回以下值:

So say there are 5 segments along a horizontal line, delimited by [0, 1.1, 2.2, 2.3, 2.8, 4] (i.e. the line goes from 0 to 4). Assume each of the segments takes an arbitrary shading values, for example, we could have the 5 shading values [0,0.88,0.55,0.11,0.44] - where 0 is black and 1 is white. Then if I wanted to plot this using 4 pixels, I would need to create 4 values, from 0...1, 1...2, etc, and would expect the calculation to return the following values for each:

0 ... 1 = 0(第一行覆盖0-> 1.1)

0...1 = 0 (this is covered by the first line segment, 0->1.1)

1 ... 2 = 0.1 * 0 + 0.9 * 0.88(第一个线段覆盖了1 ... 1.1，第二个线段覆盖了其余部分)

1...2 = 0.1 * 0 + 0.9 * 0.88 (1 ... 1.1 is covered by the first line segment, the rest by the second)

2 ... 3 = 0.2 * 0.88，0.1 * 0.55 + 0.5 * 0.11 + 0.2 * 0.44(由第二到第五条线段覆盖)

2...3 = 0.2 * 0.88, 0.1 * 0.55 + 0.5 * 0.11 + 0.2 * 0.44 (this is covered by the second to the fifth line segments )

3 ... 4 = 0.44(由最后一个线段2.8-> 4覆盖)

3...4 = 0.44 (this is covered by the last line segment, 2.8->4)

如果我想将此数据拟合为2像素长的线，则2像素将具有以下值:

Whereas if I wanted to fit this data into a 2-pixel-long line, the 2 pixels would have the following values:

0 ... 2 = 1.1/2 * 0 + 0.9/2 * 0.88

0...2 = 1.1 / 2 * 0 + 0.9 / 2 * 0.88

2 ... 4 = 0.2/2 * 0.88 + 0.1/2 * 0.55 + 0.5/2 * 0.11 + 1.2 * 0.44

2...4 = 0.2 / 2 * 0.88 + 0.1 / 2 * 0.55 + 0.5 / 2 * 0.11 + 1.2 * 0.44

这似乎是沿1d线进行下采样的正确"方法.我正在寻找一种快速实现(最好是内置的)的功能，当我沿线有(例如)一百万个点，并且只有1000个(或大约)像素可以容纳它们时.

This seems like the "right" way to do downsampling along a 1d line. I'm looking for a fast implementation (ideally something build-in) for when I have (say) a million points along the line, and only 1000 (or so) pixels to fit them in.

推荐答案

有一个纯粹的numpy解决方案.诀窍是巧妙地混合 np.searchsorted ，这会将您的常规网格放置在原始网格的最近位置，并

As you would expect, there is a purely numpy solution to this problem. The trick is to cleverly mix np.searchsorted, which will place your regular grid on the nearest bin of the original, and np.add.reduceat to compute the sums of the bins:

import numpy as np

def distribute(x, y, n):
    """
    Down-samples/interpolates the y-values of each segment across a
    domain with `n` points. `x` represents segment endpoints, so should
    have one more element than `y`. 
    """
    y = np.asanyarray(y)
    x = np.asanyarray(x)

    new_x = np.linspace(x[0], x[-1], n + 1)
    # Find the insertion indices
    locs = np.searchsorted(x, new_x)[1:]
    # create a matrix of indices
    indices = np.zeros(2 * n, dtype=np.int)
    # Fill it in
    dloc = locs[:-1] - 1
    indices[2::2] = dloc
    indices[1::2] = locs

    # This is the sum of every original segment a new segment touches
    weighted = np.append(y * np.diff(x), 0)
    sums = np.add.reduceat(weighted, indices)[::2]

    # Now subtract the adjusted portions from the right end of the sums
    sums[:-1] -= (x[dloc + 1] - new_x[1:-1]) * y[dloc]
    # Now do the same for the left of each interval
    sums[1:] -= (new_x[1:-1] - x[dloc]) * y[dloc]

    return new_x, sums / np.diff(new_x)


seg = [0, 1.1, 2.2, 2.3, 2.8, 4]
color = [0, 0.88, 0.55, 0.11, 0.44]

seg, color = distribute(seg, color, 4)
print(seg, color)

结果是

[0. 1. 2. 3. 4.] [0.    0.792 0.374 0.44 ]

这正是您在手动计算中所期望的.

Which is exactly what you expected in your manual computation.

基准

我运行了以下基准测试，以确保 EE_的解决方案和我的回答都相同，并且查看时间.我对其他解决方案进行了一些修改，使其具有与我相同的界面:

I ran the following set of benchmarks to make sure that both EE_'s solution and mine agreed on the answer, and to check out the timings. I modified the other solution slightly to have the same interface as mine:

from scipy.interpolate import interp1d

def EE_(x, y, n):
    I = np.zeros_like(x)
    I[1:] = np.cumsum(np.diff(x) * y)
    f = interp1d(x, I, bounds_error=False, fill_value=(0, I[-1]))
    pix_x = np.linspace(x[0], x[-1], n + 1)
    pix_y = (f(pix_x[1:]) - f(pix_x[:-1])) / (pix_x[1:] - pix_x[:-1])
    return pix_x, pix_y

这是测试台(方法MadPhysicist仅从上面的distribute函数重命名).对于x，输入始终为1001个元素，对于y，输入始终为1000个元素.输出数字为5、10、100、1000、10000:

And here is the testbench (the method MadPhysicist is just renamed from the distribute function above). The input is always 1001 elements for x and 1000 elements for y. The output numbers are 5, 10, 100, 1000, 10000:

np.random.seed(0x1234ABCD)

x = np.cumsum(np.random.gamma(3.0, 0.2, size=1001))
y = np.random.uniform(0.0, 1.0, size=1000)

tests = (
    MadPhysicist,
    EE_,
)

for n in (5, 10, 100, 1000, 10000):
    print(f'N = {n}')
    results = {test.__name__: test(x, y, n) for test in tests}

    for name, (x_out, y_out) in results.items():
        print(f'{name}:\n\tx = {x_out}\n\ty = {y_out}')

    allsame = np.array([[np.allclose(x1, x2) and np.allclose(y1, y2)
                         for x2, y2 in results.values()]
                        for x1, y1 in results.values()])
    print()
    print(f'Result Match:\n{allsame}')

    from IPython import get_ipython
    magic = get_ipython().magic

    for test in tests:
        print(f'{test.__name__}({n}):\n\t', end='')
        magic(f'timeit {test.__name__}(x, y, n)')

我将跳过数据和协议打印输出(结果相同)，并显示时间:

I will skip the data and agreement printouts (the results are identical), and show the timings:

N = 5
MadPhysicist: 50.6 µs ± 349 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
EE_:           110 µs ± 568 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

N = 10
MadPhysicist: 50.5 µs ± 732 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
EE_:           111 µs ± 635 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)   

N = 100
MadPhysicist: 54.5 µs ± 284 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
EE_:           114 µs ± 215 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

N = 1000
MadPhysicist: 107 µs ± 5.73 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
EE_:          148 µs ± 5.11 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

N = 10000
MadPhysicist: 458 µs ± 2.21 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
EE_:          301 µs ± 4.57 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

您可以看到，在较小的输出尺寸下，numpy解决方案要快得多，这可能是因为开销占主导.但是，在更多的断点处，scipy解决方案变得更快.您必须比较不同的输入大小，才能真正了解时序的计算原理，而不仅仅是不同的输出大小.

You can see that at smaller output sizes, the numpy solution is much faster, probably because the overhead dominates. However, at larger numbers of break points, the scipy solution becomes much faster. You would have to compare at different input sizes to get a real idea of how the timing works out, not just different output sizes.

这篇关于numpy:大量线段/点的快速有规律的平均间隔的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！