Эффективно применить функцию к сферической окрестности в массиве NumPy - PullRequest
Новая
       13

Эффективно применить функцию к сферической окрестности в массиве NumPy

0 голосов
/ 04 декабря 2018

У меня есть трехмерный массив значений с плавающей точкой в ​​Python.Мне нужно получить все элементы в сфере радиуса r, начиная с центральной точки P (x, y, z).Затем я хочу применить к точкам сферы функцию, которая обновляет их значения и для этого требуется расстояние до центральной точки.Я делаю эти шаги много раз и для больших значений радиуса, поэтому я хотел бы найти решение, которое было бы максимально эффективным.

Мое текущее решение проверяет только точки в ограничительной рамке сферы, как указано здесь: Использование QuadTree для получения всех точек в пределах ограничительной окружности .Эскиз кода выглядит следующим образом: 

# P(x, y, z): center of the sphere
for k1 in range(x - r, x + r + 1):
    for k2 in range(y - r, y + r + 1):
        for k3 in range(z - r, z + r + 1):
            # Sphere center - current point distance
            dist = np.sum((np.array([k1, k2, k3]) - np.array([x, y, z])) ** 2)

            if (dist <= r * r):
                # computeUpdatedValue(distance, radius): function that computes the new value of the matrix in the current point 
                newValue = computeUpdatedValue(dist, r)

                # Update the matrix
                mat[k1, k2, k3] = newValue

Однако я подумал, что применение маски для извлечения точек, а затем их обновления на основе расстояния в векторизованном виде более эффективно.Я видел, как применить кольцевое ядро ​​( Как применить маску в форме диска к массиву numpy? ), но я не знаю, как эффективно применять функцию (в зависимости от индексов) к каждому изэлементы маски.

1 Ответ

0 голосов
/ 04 декабря 2018

РЕДАКТИРОВАТЬ: Если ваш массив очень большой по сравнению с регионом, который вы обновляете, решение ниже займет гораздо больше памяти, чем необходимо.Вы можете применить ту же идею, но только к области, где сфера может упасть: 

def updateSphereBetter(mat, center, radius):
    # Find beginning and end of region of interest
    center = np.asarray(center)
    start = np.minimum(np.maximum(center - radius, 0), mat.shape)
    end = np.minimum(np.maximum(center + radius + 1, 0), mat.shape)
    # Slice region of interest
    mat_sub = mat[tuple(slice(s, e) for s, e in zip(start, end))]
    # Center coordinates relative to the region of interest
    center_rel = center - start
    # Same as before but with mat_sub and center_rel
    ind = np.indices(mat_sub.shape)
    ind = np.moveaxis(ind, 0, -1)
    dist_squared = np.sum(np.square(ind - center_rel), axis=-1)
    mask = dist_squared <= radius * radius
    mat_sub[mask] = computeUpdatedValue(dist_squared[mask], radius)

Обратите внимание, что, поскольку mat_sub является представлением mat, при обновлении он обновляет исходный массив, поэтомутот же результат, что и раньше, но с меньшими ресурсами.

Вот небольшое доказательство концепции.Я определил computeUpdatedValue так, чтобы он показывал расстояние от центра, а затем построил несколько «отрезков» примера: 

import numpy as np
import matplotlib.pyplot as plt

def updateSphere(mat, center, radius):
    # Make array of all index coordinates
    ind = np.indices(mat.shape)
    # Compute the squared distances to each point
    ind = np.moveaxis(ind, 0, -1)
    dist_squared = np.sum(np.square(ind - center), axis=-1)
    # Make a mask for squared distances within squared radius
    mask = dist_squared <= radius * radius
    # Update masked values
    mat[mask] = computeUpdatedValue(dist_squared[mask], radius)

def computeUpdatedValue(dist_squared, radius):
    # 1 at the center of the sphere and 0 at the surface
    return np.clip(1 - np.sqrt(dist_squared) / radius, 0, 1)

mat = np.zeros((100, 60, 80))
updateSphere(mat, [50, 20, 40], 20)

plt.subplot(131)
plt.imshow(mat[:, :, 30], vmin=0, vmax=1)
plt.subplot(132)
plt.imshow(mat[:, :, 40], vmin=0, vmax=1)
plt.subplot(133)
plt.imshow(mat[:, :, 55], vmin=0, vmax=1)

Вывод:

Sphere update

...