title: 3D Median box filter on the 4-th dimension of a point cloud (Python)
body: <p>I'm looking for an efficient way to compute a moving 3D median filter of size m<em>n</em>p on <strong>the 4th dimension</strong> "U" of a large 3D point cloud (x,y,z). This differs from typical filter for Nd images in two ways: (1) there can be several values in each image voxel (several points falling in that voxels) and (2) Some voxels can be empty, and should be replaced by the median in the box around them.</p>
  
  <p>I have converged to the following method, which however seems not so efficient when the number of points increases (while it performs better for increasing box averaging windows) ...</p>
  
  <pre><code>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&gt;=0)&amp;(neigh&lt;=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]&gt;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
  </code></pre>
  
  <p>The time statistics are the following :</p>
  
  <pre><code>filt_size       10^5 points     10^6 points     10^7 points
  [1,1,1]         10.89 s         40.09 s         137.97
  [2,1,1]         13.64 s         68.46 s         x
  [3,1,1]         16.67 s         94.10 s         x
  [2,2,2]         25.57 s         207.75 s
  </code></pre>
  
  <p>Some ideas to improve the methods ?
  Thank you very much</p>