Суммируйте ndarray по 2d массиву в Python - PullRequest
Новая
       18

Суммируйте ndarray по 2d массиву в Python

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

Я хочу суммировать трехмерный массив dat, используя индексы, содержащиеся в двумерном массиве idx.

Рассмотрим пример ниже. Для каждого поля по dat[:, :, i] я хочу вычислить медиану по некоторому индексу idx. Желаемый вывод (out) - это двумерный массив, строки которого записывают индекс, а столбцы - поле. Следующий код работает, но не очень эффективен. Есть предложения? 

import numpy as np
dat = np.arange(12).reshape(2, 2, 3)
idx = np.array([[0, 0], [1, 2]])

out = np.empty((3, 3))
for i in np.unique(idx):
    out[i,] = np.median(dat[idx==i], axis = 0)
print(out)

Выход: 

[[ 1.5  2.5  3.5]
 [ 6.   7.   8. ]
 [ 9.  10.  11. ]]

1 Ответ

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

Чтобы лучше представить проблему, я буду ссылаться на размеры массива 2x2 как строки и столбцы, а на размерность 3 - на глубину. Я буду ссылаться на векторы вдоль третьего измерения как на «пиксели» (пиксели имеют длину 3), а плоскости вдоль первых двух измерений - на «каналы».

Ваш цикл накапливает набор пикселей, выбранных маской idx == i, и получает медиану каждого канала в этом наборе. В результате получается массив Nx3, где N - это количество различных приливов, которые у вас есть.

Однажды, обобщенные уфунки будут вездесущими в numpy, и np.median будет такой функцией. В этот день вы сможете использовать reduceat magic 1 , чтобы сделать что-то вроде 

unq, ind = np.unique(idx, return_inverse=True)
np.median.reduceat(dat.reshape(-1, dat.shape[-1]), np.r_[0, np.where(np.diff(unq[ind]))[0]+1])

1 См. Применение операции к неравномерно разделенным частям массива для получения дополнительной информации о конкретном типе магии.

Поскольку в настоящее время это невозможно, вместо этого можно использовать scipy.ndimage.median. Эта версия позволяет вычислять медианы по набору отмеченных областей в массиве, что в точности соответствует idx. Этот метод предполагает, что ваш индексный массив содержит N плотно упакованных значений, все из которых находятся в range(N). В противном случае операции изменения формы не будут работать должным образом.

Если это не так, начните с преобразования idx: 

_, ind = np.unique(idx, return_inverse=True)
idx = ind.reshape(idx.shape)

OR 

idx = np.unique(idx, return_inverse=True)[1].reshape(idx.shape)

Поскольку вы фактически вычисляете отдельную медиану для каждого региона и канала, вам потребуется набор меток для каждого канала. Уточните idx, чтобы иметь отдельный набор индексов для каждого канала: 

chan = dat.shape[-1]
offset = idx.max() + 1
index = np.stack([idx + i * offset for i in range(chan)], axis=-1)

Теперь index имеет идентичный набор областей, определенных в каждом канале, который вы можете использовать в scipy.ndimage.median: 

out = scipy.ndimage.median(dat, index, index=range(offset * chan)).reshape(chan, offset).T

Входные метки должны быть плотно упакованы от нуля до offset * chan, чтобы index=range(offset * chan) работал правильно, а операция reshape - для правильного количества элементов. Окончательная транспонирование - просто артефакт того, как расположены метки.

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

import numpy as np
from scipy.ndimage import median

dat = np.arange(12).reshape(2, 2, 3)
idx = np.array([[0, 0], [1, 2]])

def summarize(dat, idx):
    idx = np.unique(idx, return_inverse=True)[1].reshape(idx.shape)
    chan = dat.shape[-1]
    offset = idx.max() + 1
    index = np.stack([idx + i * offset for i in range(chan)], axis=-1)
    return median(dat, index, index=range(offset * chan)).reshape(chan, offset).T

print(summarize(dat, idx))
...