python-将函数有效地应用于numpy数组中的球面邻域
作者:互联网
我在Python中有一个浮点值的3D numpy数组.
我需要检索半径为r的球体中的所有元素,从
中心点P(x,y,z).然后,我要在球体上应用点函数
更新其值,并需要到中心点的距离才能执行此操作.我经常执行这些步骤,并且
大半径值,所以我想有一个同样有效的解决方案
尽可能.
我当前的解决方案仅检查球体边界框中的点,
如此处所示:Using a QuadTree to get all points within a bounding circle.
代码草图如下所示:
# 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
但是,我认为应用遮罩来检索点,然后
基于距离以矢量化方式更新它们的效率更高.
我已经看到了如何应用循环内核
(How to apply a disc shaped mask to a numpy array?),
但是我不知道如何在蒙版的每个元素上有效地应用该函数(取决于索引).
解决方法:
编辑:如果与正在更新的区域相比,您的数组很大,则下面的解决方案将占用比必要更多的内存.您可以将相同的想法应用到球体可能掉落的区域:
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)
输出:
标签:vectorization,mask,python,numpy 来源: https://codeday.me/bug/20191108/2007544.html