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

Question

我有一个二维数据数组，我试图以有效的方式获得有关其中心值的轮廓。所以输出应该是两个一维数组：一个是距离中心的距离值，另一个是原始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需要等距箱，但情况并非如此。任何建议？

Answer 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 

Answer 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 

但我还是不愉快的事实，我不得不使用列表理解（或蟒蛇循环），因为它对于非常大的矩阵可能会很慢。

Answer 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张图片上交叉。