我正在寻找一种有效的方法来计算一个移动的三维中值滤波器，其大小为mnp on一个大型三维点云(x，y，z)的第四维“U”。这与Nd图像的典型过滤器有两个不同之处：(1)每个图像体素中可能有多个值(多个点落在该体素中)和(2)某些体素可以为空，并应替换为其周围框中的中值。在
我已经收敛到下面的方法，但是当点数增加时，这种方法似乎不那么有效(虽然它在增加框平均窗口时表现更好)。。。在import itertools
import numpy as np
import time
import sys
npts=10**5 # Number of points in the cloud
L=100 # Dimension of the domain (grid)
filt_size=[1,1,1] # extent of the median filter on each size of the grid voxel
# Generate random points and scalar u
x=np.random.random(npts)*L
y=np.random.random(npts)*L
z=np.random.random(npts)*L
u=np.random.random(npts)
# Initiate time clock
t0=time.clock()
# Take grid position of points
x=np.uint32(np.round(x))
y=np.uint32(np.round(y))
z=np.uint32(np.round(z))
# Index point in a flatten way
index=x*L*L+y*L+z
idmax=np.max(index) # Maximum for latter use
idsort=np.argsort(index) # Sort indexes to be quicker ?
index=index[idsort]
u=np.float16(u[idsort])
# Group by grid pixel (unique indexes given by udi and their extent by ui)
idu,ui=np.unique(index,return_index=True)
# Neighboor indexes of the index 0
neighboors=np.array([i*L*L+j*L+k for i, j,k in itertools.product(range(-filt_size[0],filt_size[0]+1), range(-filt_size[1],filt_size[1]+1),range(-filt_size[2],filt_size[2]+1))],dtype='int32')
# Make L^3 list of empty arrays in which we will aggregate all points in the neighboorood of a voxel
U=[np.array([],dtype='float16')]*int((L+1)*(L+1)*(L+1))
for i in range(len(idu)-1):
ul=u[ui[i]:ui[i+1]] # the scalar data of voxel idu[i]
neigh=neighboors+idu[i] # Neighbooring cells TO where the scalar data should be copied
for k in neigh[(neigh>=0)&(neigh<=idmax)]: # Take only neighboors in the grid domain
U[k]=np.concatenate((ul,U[k])) # stack with existent data
U[k].sort(axis=0) # sort
sys.stdout.write('\r progress {:d}%'.format(i*100/len(idu)))
sys.stdout.flush()
# Finally, gather the median by taking the middle element of each array
# empty voxels are filled by nan
Umed=np.zeros((L+1)**3,dtype='float16')+np.nan
for i in range(len(U)):
if U[i].shape[0]>0: # Check if not-empty
Umed[i]=U[i][np.round(U[i].shape[0]/2)]
sys.stdout.write('\r progress {:d}%'.format(i*100/len(U)))
sys.stdout.flush()
# Reshape Umed in a (L+1)*(L+1)*(L+1) array
Umed=Umed.reshape((L+1,L+1,L+1))
print 'Elapsed :',time.clock()-t0
时间统计如下：
^{pr2}$
改进方法的一些想法？
非常感谢