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

Question

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

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

考虑这些矩阵

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，这非常方便(但对我不起作用)。

我一直在尝试将其扩展为 float 距离，就像我的版本 r 一样，但到目前为止还没有成功。 bincount 不适用于 float ，histogram 需要等间距的 bin，但事实并非如此。有什么建议吗？

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