具有浮點索引的二維矩陣的徑向輪廓

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張圖片上交叉。