具有浮点索引的二维矩阵的径向轮廓答案

【问题标题】：Radial profile of 2D matrix with float indexes具有浮点索引的二维矩阵的径向轮廓
【发布时间】：2017-07-28 07:46:05
【问题描述】：

我有一个 2D 数据数组，我正在尝试以有效的方式获取关于其中心的值的配置文件。所以输出应该是两个一维数组：一个是到中心的距离值，另一个是原始 2D 中距离中心该距离的所有值的平均值。

每个索引与中心的距离都不是整数，这使我无法使用一些已知的解决方案来解决该问题。请允许我解释一下。

考虑这些矩阵

data = np.random.randn(5,5)
L = 2
x = np.arange(-L,L+1,1)*2.5
y = np.arange(-L,L+1,1)*2.5
xx, yy = np.meshgrid(x, y)
r = np.sqrt(xx**2. + yy**2.)

所以矩阵是

In [30]: r
Out[30]: 
array([[ 7.07106781,  5.59016994,  5.        ,  5.59016994,  7.07106781],
       [ 5.59016994,  3.53553391,  2.5       ,  3.53553391,  5.59016994],
       [ 5.        ,  2.5       ,  0.        ,  2.5       ,  5.        ],
       [ 5.59016994,  3.53553391,  2.5       ,  3.53553391,  5.59016994],
       [ 7.07106781,  5.59016994,  5.        ,  5.59016994,  7.07106781]])

In [31]: data
Out[31]: 
array([[ 1.27603322,  1.33635284,  1.93093228,  0.76229675, -0.00956535],
       [ 0.69556071, -1.70829753,  1.19615919, -1.32868665,  0.29679494],
       [ 0.13097791, -1.33302719,  1.48226442, -0.76672223, -1.01836614],
       [ 0.51334771, -0.83863115, -0.41541794,  0.34743342,  0.1199237 ],
       [-1.02042539,  0.90739383, -2.4858624 , -0.07417987,  0.90748933]])

在这种情况下，距离索引的预期输出应该是 array([ 0. , 2.5 , 3.53553391, 5. , 5.59016994, 7.07106781])，第二个长度相同的数组应该是那些对应距离处所有值的平均值：array([ 0.98791323, -0.32496927, 0.37221219, -0.6209728 , 0.27986926, 0.04060628])。

来自this answer 有一个非常好的函数可以计算任意点的轮廓。然而，他的方法的问题在于它通过索引距离来近似距离r。所以他的r 对我来说是这样的：

array([[2, 2, 2, 2, 2],
       [2, 1, 1, 1, 2],
       [2, 1, 0, 1, 2],
       [2, 1, 1, 1, 2],
       [2, 2, 2, 2, 2]])

这对我来说是一个很大的不同，因为我正在使用小矩阵。然而，这个近似值允许他使用np.bincount，这非常方便（但对我不起作用）。

我一直在尝试将其扩展为浮动距离，例如我的版本r，但到目前为止还没有运气。 bincount 不适用于浮点数，histogram 需要等间距的 bin，但事实并非如此。有什么建议吗？

【问题讨论】：

使用((xx**2. + yy**2.)/6.25).astype(int) 作为bincount 的bin 怎么样？
或者在r上使用np.digitize。
@PaulPanzer 我不明白如何使用数字化来做到这一点。愿意举个例子吗？
这类似于@Divakar 的库函数所做的。您将提供距离的升序序列来定义要组合在一起的环。将这些 bin 应用于r，您将得到一个与 Divakar 的 Out[280] 非常相似的矩阵。 -- 不相关：我想知道您是否要将相同的距离配置文件应用于许多图像？
@PaulPanzer 给定集合的距离配置文件是固定的，所以我只需将其与图像列表一起传递给函数

标签： python performance numpy matrix

【解决方案1】：

方法#1

def radial_profile_app1(data, r):
    mid = data.shape[0]//2
    ids = np.rint((r**2)/r[mid-1,mid]**2).astype(int).ravel()
    count = np.bincount(ids)

    R = data.shape[0]//2 # Radial profile radius
    R0 = R+1
    dists = np.unique(r[:R0,:R0][np.tril(np.ones((R0,R0),dtype=bool))])

    mean_data = (np.bincount(ids, data.ravel())/count)[count!=0]
    return dists, mean_data

对于给定的样本数据-

In [475]: radial_profile_app1(data, r)
Out[475]: 
(array([ 0.        ,  2.5       ,  3.53553391,  5.        ,  5.59016994,
         7.07106781]),
 array([ 1.48226442  , -0.3297520425, -0.8820454775, -0.3605795875,
         0.5696863263,  0.2883829525]))

方法 #2

def radial_profile_app2(data, r):
    R = data.shape[0]//2 # Radial profile radius
    range_arr = np.arange(-R,R+1)
    ids = (range_arr[:,None]**2 + range_arr**2).ravel()
    count = np.bincount(ids)

    R0 = R+1
    dists = np.unique(r[:R0,:R0][np.tril(np.ones((R0,R0),dtype=bool))])

    mean_data = (np.bincount(ids, data.ravel())/count)[count!=0]
    return dists, mean_data

运行时测试 -

In [562]: # Setup inputs
     ...: N = 2001
     ...: data = np.random.randn(N,N)
     ...: L = (N-1)//2
     ...: x = np.arange(-L,L+1,1)*2.5
     ...: y = np.arange(-L,L+1,1)*2.5
     ...: xx, yy = np.meshgrid(x, y)
     ...: r = np.sqrt(xx**2. + yy**2.)
     ...: 

In [563]: out01, out02 = radial_profile_app1(data, r)
     ...: out11, out12 = radial_profile_app2(data, r)
     ...: 
     ...: print np.allclose(out01, out11)
     ...: print np.allclose(out02, out12)
     ...: 
True
True

In [566]: %timeit radial_profile_app1(data, r)
     ...: %timeit radial_profile_app2(data, r)
     ...: 
10 loops, best of 3: 114 ms per loop
10 loops, best of 3: 91.2 ms per loop

【讨论】：

看起来很有趣，但我没有看到一种简单的方法将所有内容都作为原始 r 函数的函数。
@TomCho 那么Approach #3 可能会满足您的需求。
我想了解你现在做了什么。我想你可能误解了我在找什么。就我而言，我正在寻找的输出具有array([ 0. , 2.5 , 3.53553391, 5. , 5.59016994, 7.07106781])，并且值是means 上具有这些索引的data 矩阵上的每个值。如果这不是您所理解的，请告诉我，以便我更新我的答案，
@TomCho 是的，我完全错过了在问题或 cmets 的任何地方看到任何提及 mean 或 average 一词。所以，请随时更新相同的问题:)
我刚刚发布了一个编辑，希望它更容易理解。但是您的回答给了我一个可能解决我的问题的想法，尽管它并不是非常优化。如果它有效，我会在一秒钟内发布它

【解决方案2】：

得到了我对这个功能的期望：

def radial_prof(data, r):
    uniq = np.unique(r)
    prof = np.array([ np.mean(data[ r==un ]) for un in uniq ])
    return uniq, prof

但我仍然对必须使用列表解析（或 python 循环）这一事实不满意，因为对于非常大的矩阵来说它可能会很慢。

【讨论】：

【解决方案3】：

这是一种间接排序方法，如果批量大小和/或箱数很大，它应该可以很好地扩展。排序是 O(n log n) 所有直方图是 O(n)。我还添加了一些不科学的速度测试。对于速度测试，我使用平面索引，但我保留了 2d 索引代码，因为它在处理不同尺寸的图像等时更灵活。

import numpy as np

# this need only be run once per batch
def r_to_ind(r, dist_bins="auto"):
    f = np.argsort(r.ravel())
    if dist_bins == "auto":
        rs = r.ravel()[f]
        bins = np.where(np.r_[True, rs[1:]!=rs[:-1]])[0]
        dist_bins = rs[bins]
    else:
        bins = np.searchsorted(r.ravel()[f], dist_bins)
    denom = np.diff(np.r_[bins, r.size])
    return f, np.unravel_index(f, r.shape), bins, denom, dist_bins

# this is with adjustable offset
def profile_xy(image, yx, ij, bins, nynx, denom):
    (y, x), (i, j), (ny, nx) = yx, ij, nynx
    return np.add.reduceat(image[i + y - ny//2, j + x - nx//2], bins) / denom

# this is fixed
def profile_xy_no_offset(image, ij, bins, denom):
    return np.add.reduceat(image[ij], bins) / denom

# this is fixed and flat
def profile_xy_no_offset_flat(image, k, bins, denom):
    return np.add.reduceat(image.ravel()[k], bins) / denom

data = np.array([[ 1.27603322,  1.33635284,  1.93093228,  0.76229675, -0.00956535],
       [ 0.69556071, -1.70829753,  1.19615919, -1.32868665,  0.29679494],
       [ 0.13097791, -1.33302719,  1.48226442, -0.76672223, -1.01836614],
       [ 0.51334771, -0.83863115, -0.41541794,  0.34743342,  0.1199237 ],
       [-1.02042539,  0.90739383, -2.4858624 , -0.07417987,  0.90748933]])

r = np.array([[ 7.07106781,  5.59016994,  5.        ,  5.59016994,  7.07106781],
       [ 5.59016994,  3.53553391,  2.5       ,  3.53553391,  5.59016994],
       [ 5.        ,  2.5       ,  0.        ,  2.5       ,  5.        ],
       [ 5.59016994,  3.53553391,  2.5       ,  3.53553391,  5.59016994],
       [ 7.07106781,  5.59016994,  5.        ,  5.59016994,  7.07106781]])

f, (i, j), bins, denom, dist_bins = r_to_ind(r)

result = profile_xy(data, (2, 2), (i, j), bins, (5, 5), denom)
print(dist_bins)
# [ 0.          2.5         3.53553391  5.          5.59016994  7.07106781]
print(result)
# [ 1.48226442 -0.32975204 -0.88204548 -0.36057959  0.56968633  0.28838295]

#########################

from timeit import timeit

n = 2001
batch = 100
fake = 10

a = np.random.random((fake, n, n))
l = np.linspace(-1, 1, n)**2
r = sum(np.ix_(l, l))

def run_all():
    f, ij, bins, denom, dist_bins = r_to_ind(r)
    for b in range(batch):
        profile_xy_no_offset_flat(a[b%fake], f, bins, denom)

print(timeit(run_all, number=10))
# 47.4157 (for 10 batches of 100 images of size 2001x2001)
# and my computer is slower than Divakar's ;-)

我做了更多的基准测试，将我的方法与@Divakar 的方法 3 进行比较，将所有可预计算的东西都剥离成每批运行一次的函数。一般发现：它们与我的相似，前期成本更高，但速度更快。但它们每批只能交叉 100 张左右的图片。

【讨论】：