Вычисление функции эффективной плотности

Question

У меня есть большое изображение в форме массива numpy (opencv возвращает его как 2d-массив из 3 значений uint8) и хочет вычислить сумму гауссовых ядер для каждого пикселя, то есть (в LaTeX поддержка по-прежнему отсутствует SO есть?):

для N различных ядер с заданной массой w, средней и диагональной ковариационной матрицей.

Так что я хочу функцию compute_densities(image, kernels) -> numpy array of floats. Каков наилучший способ сделать это эффективно в python? Я был бы удивлен, если бы не было библиотечной функции в scipy для этого, но у меня была статистика в uni давным-давно, поэтому я немного запутался в деталях документации.

В принципе Я хочу следующее, просто более эффективное, чем наивный питон (2pi^{- 3/2} игнорируется, поскольку это постоянный фактор, который не имеет для меня значения, поскольку меня интересуют только отношения между вероятностями)

def compute_probabilities(img, kernels): 
    np.seterr(divide='ignore') # 1/covariance logs an error otherwise 
    result = np.zeros((img.shape[0], img.shape[1])) 
    for row_pos, row_val in enumerate(img): 
     for col_pos, val in enumerate(row_val): 
      prob = 0.0 
      for kernel in kernels: 
       mean, covariance, weight = kernel 
       val_sub_mu = np.array([val]).T - mean 
       cov_inv = np.where(covariance != 0, 1/covariance, 0) 
       tmp = val_sub_mu.T.dot(cov_inv).dot(val_sub_mu) 
       prob += weight/np.sqrt(np.linalg.norm(covariance)) * \ 
         math.exp(-0.5 * tmp) 
      result[row_pos][col_pos] = prob 
    np.seterr(divide='warn') 
    return result 

Вход: cv2.imread на некоторый jpg, который дает 2d массив (высота x ширина) структуры 3 uint8, содержащей 3 цветовых канала.

Ядра является namedtuple('Kernel', 'mean covariance weight'), значит, вектор, ковариация является 3x3 матрица со всем, но диагональ равна нулю и вес поплавка 0 < weight < 1. Для простоты я только указать диагонали, а затем преобразовать его в матрицу 3х3 впоследствии: (представление не высечены на камне не все равно, как это представлено так быть свободным, чтобы изменить все, что):

some_kernels = [ 
    Kernel(np.array([(73.53, 29.94, 17.76)]), np.array([(765.40, 121.44, 112.80)]), 0.0294), 
    ... 
] 

def fixup_kernels(kernels): 
    new_kernels = [] 
    for kernel in kernels: 
     cov = np.zeros((3, 3)) 
     for pos, c in enumerate(kernel.covariance[0]): 
      cov[pos][pos] = c 
     new_kernels.append(Kernel(kernel.mean.T, cov, kernel.weight)) 
    return new_kernels 

some_kernels = fixup_kernels(some_kernels) 
img = cv2.imread("something.jpg") 
result = compute_probabalities(img, some_kernels) 

Answer 1

Как я понимаю, ваша цель - оценить плотность вероятности каждого пикселя изображения по отношению к смеси многомерных нормальных распределений.

[Существует в настоящее время (2013-11-20) ошибка в коде в вашем вопросе и ответе @Alex I-х - | | окружающих \Sigma в приведенном выше уравнении фактически означает определитель, а не нормы вектора - см например here. В случае диагональной ковариации детерминант является просто произведением диагональных элементов.]

Очень удобно реализовать расчет плотности с точки зрения операций с массивом numpy. Следующая реализация использует сферические (то есть по диагонали) природы матриц ковариации в вашей задаче:

def compute_probabilities_faster(img, kernels): 
    means, covs, weights = map(np.dstack, zip(*kernels)) 
    pixels_as_rows = img.reshape((-1, 3, 1)) 
    responses = np.exp(-0.5 * ((pixels_as_rows - means) ** 2/covs).sum(axis=1)) 
    factors = 1./np.sqrt(covs.prod(axis=1) * ((2 * np.pi) ** 3)) 
    return np.sum(responses * factors * weights, axis=2).reshape(img.shape[:2]) 

Эта функция действует непосредственно на ядрах, как они первоначально представлены , то есть. без модификация fixup_kernels функция. Когда коэффициент нормализации (2 * np.pi) ** 3 удаляется (и звонки на linalg.norm заменяются на linalg.det), эта функция соответствует выходу вашего кода (достаточно для удовлетворения np.allclose).

Ближайшая функциональность вне коробки в SciPy (как 0,13) является реализацией ядра оценивания плотности в scipy.stats (см here), который определяет очень похожий распределения где матрицы ковариаций являются одинаковыми для каждое ядро ​​- по этой причине не подходит для вашей проблемы.

Answer 2

Способ получения производительности от Python - не использовать Python.

Существует много пакетов, которые работают с синтаксисом Python, но затем используют C или C++. Сам NumPy делает это. Ваша проблема, похоже, специально разработана для чего-то вроде Cython или numexpr. Обе эти ссылки показывают вам, как использовать любую систему для ядер на векторах NumPy.

Редактировать: Я бы хотел, чтобы один из моих подставных лиц дал мне знать, как я ошибаюсь. Я предлагаю вам пойти, если не удалось найти предварительно созданную функцию. Если вы знаете способ сделать это с большей производительностью, чем что-то вроде Cython или numexpr (т. Е. Способ писать C в синтаксисе Python), я бы хотел его услышать.

Answer 3

EDIT

Я проверил, что это дает те же результаты, что и исходный код:

def compute_probabilities_fast(img, kernels): 
    np.seterr(divide='ignore') 
    result = np.zeros((img.shape[0], img.shape[1])) 
    for kernel in kernels: 
     mean, covariance, weight = kernel 
     cov_inv = np.where(covariance != 0, 1/covariance, 0) 
     mean = mean[:,0] 
     img_sub_mu = img - mean 
     img_tmp = np.sum(img_sub_mu.dot(cov_inv) * img_sub_mu, axis=2) 
     result += (weight/np.sqrt(np.linalg.norm(covariance))) * np.exp(-0.5 * img_tmp) 
    return result 

Объяснение:

mean[:,0] делает форму просто (3,) вместо (3, 1).

img - mean транслирует на весь образ и вычитает средства из каждого пикселя.

img_sub_mu.dot(cov_inv) приблизительно эквивалентен val_sub_mu.T.dot(cov_inv).

np.sum(... * img_sub_mu, axis=2) приблизительно эквивалентен .dot(val_sub_mu). Нельзя использовать точку, хотя, потому что это добавит дополнительные размеры.Например, массив M x N x K, усеянный массивом M x K x N, приведет к результату M x N x M x N, точка ведет себя по-разному на одномерные и многомерные данные. Таким образом, мы просто умножаем элемент и затем суммируем по последнему измерению.

На самом деле «сумма гауссовых ядер» в этом вопросе меня смутила. Запрошенным является вычисление, в котором для каждого выходного пикселя значение зависит только от только от от входного значения для того же пикселя, но не от значений соседних пикселей. Таким образом, это не что иное, как гауссовское размытие (которое будет использовать свертку), это всего лишь расчет, выполняемый по каждому пикселю в отдельности.

P.S. 1/covariance проблематично. Вы уверены, что не хотите np.linalg.inv(covariance)?

OLD ОТВЕТ

Это похоже на то, что вы хотите, один из них:

scipy.signal.convolve2d

scipy.ndimage.filters.convolve

Вопрос немного запутанным, вы пытаетесь вычислить кучу изображений, свернутых с разными гауссами, или единое изображение, свернутое суммой гауссиан? Разделяется ли ваше ядро? (Если да, используйте две свертки Mx1 и 1xN вместо одного MxN). Функции scipy, которые вы будете использовать, будут одинаковыми в любом случае.

Также, конечно, вы должны предварительно вычислить свои ядра с помощью комбинации numpy.random.normal и meshgrid.