Создание быстрого пользовательского фильтра для объемных 3D-массивов - PullRequest
0 голосов
/ 26 июня 2018

В целом, я хочу создать фильтр для вычисления среднего значения круговых величин для трехмерных массивов.

Я заглянул в scipy.ndimage.generic_filter, но не смог скомпилировать фильтр, как описано в https://ilovesymposia.com/tag/numba,, по-видимому, из-за ошибки в numba в Windows.

Затем я попытался создать свою собственную реализацию, которая проходит по массиву, и надеялся, что сможет потом его объединить. Он работает нормально (и медленно) без numba, но компиляция jit не удалась, и я не могу декодировать TypingError.

сетка numpy не поддерживается, поэтому поведение этого тоже должно было быть построено (дешевая версия).

from numba import njit
import numpy as np

@njit
def my_meshgrid(i_, j_,k_):
    #Note: axes 0 and 1 are swapped!
    shape = (len(j_), len(i_), len(k_))
    io = np.empty(shape, dtype=np.int32)
    jo = np.empty(shape, dtype=np.int32)
    ko = np.empty(shape, dtype=np.int32)
    for i in range(len(i_)):
        for j in range(len(j_)):
            for k in range(len(k_)):
                io[j,i,k] = i_[i]
                jo[j,i,k] = j_[j]
                ko[j,i,k] = k_[k]
    return [io,jo, ko]
t3 = my_meshgrid(range(5), range(5,7), range(7,10))
#

@njit
def get_footprint(arr, i , j , k, size=3):
    s = size
    ranges = [range(d-s+1+1,d+s-1) for d in [i,j,k]]
    #Mirror the case where indexes are less than zero
    ind = np.abs(np.meshgrid(*ranges))
    #Mirror the case where indexes are higher than arr.shape:
    for d in range(len(arr.shape)):
        indd = ind[d] - arr.shape[d]
        indd *= -1
        indd = np.abs(indd)
        indd *= -1
        ind[d] = indd
    return arr[ind]


@njit
def mean_angle_filter(degrees, size = 3):
    size = [size]*len(degrees.shape)
    out = np.empty_like(degrees)
    for i in range(degrees.shape[0]):
        for j in range(degrees.shape[1]):
            for k in range(degrees.shape[2]):
                out[i,j,k] = mean_angle(get_footprint(degrees, i,j,k,3))
    return out


@njit
def mean_angle(degrees):
    '''
    https://en.wikipedia.org/wiki/Mean_of_circular_quantities
    '''
    x = np.mean(np.cos(degrees*np.pi/180))
    y = np.mean(np.sin(degrees*np.pi/180))
    return np.arctan2(y,x)*180/np.pi


degrees = np.random.random([20]*3)*90
mean_angle_filter(degrees)

Будучи новичком в numba, я был бы рад получить исправление для этой (или подобной) реализации, но любая (быстрая) реализация фильтра mean_angle в numpy также была бы признательна

1 Ответ

0 голосов
/ 26 июня 2018

Вы можете значительно упростить свой код:

  1. Зеркальное отображение на границах может быть выполнено с помощью np.pad()
  2. Перекрывающиеся окна данных могутэффективно сделать это с помощью SciKit-Image skimage.util.view_as_windows().
  3. Вычисление среднего по оси можно выполнить с помощью np.mean(..., axis=x), где x - это int или tupleиз int s, обозначающих оси, которые вы хотите иметь в виду.

Соберите все это вместе, и вы получите очень простую векторизованную реализацию, такую ​​как

import numpy as np
import skimage

def mean_angle(degrees, axis=None):
    '''
    https://en.wikipedia.org/wiki/Mean_of_circular_quantities
    '''
    x = np.mean(np.cos(degrees*np.pi/180), axis=axis)
    y = np.mean(np.sin(degrees*np.pi/180), axis=axis)
    return np.arctan2(y,x)*180/np.pi


out = mean_angle(
    skimage.util.view_as_windows(
        np.pad(degrees, (1, 1), mode='symmetric'),
        window_shape=(3, 3, 3)
    ),
    axis=(-1, -2, -3)
)
...