<p>经过一些优化之后,我们为一个通用的3D图像提供了这个功能:</p>
<pre><code>def variance_filter( img, VAR_FILTER_SIZE ):
from numpy.lib.stride_tricks import as_strided
WIN_SIZE=(2*VAR_FILTER_SIZE)+1
if ~ VAR_FILTER_SIZE%2==1:
print 'Warning, VAR_FILTER_SIZE must be ODD Integer number '
# hack -- this could probably be an input to the function but Alessandro is lazy
WIN_DIMS = [ WIN_SIZE, WIN_SIZE, WIN_SIZE ]
# Check that there is a 3D image input.
if len( img.shape ) != 3:
print "\t variance_filter: Are you sure that you passed me a 3D image?"
return -1
else:
DIMS = img.shape
# Set up a windowed view on the data... this will have a border removed compared to the img_in
img_strided = as_strided(img, shape=(DIMS[0]-WIN_DIMS[0]+1, DIMS[1]-WIN_DIMS[1]+1, DIMS[2]-WIN_DIMS[2]+1, WIN_DIMS[0], WIN_DIMS[1], WIN_DIMS[2] ), strides=img.strides*2)
# Calculate variance, vectorially
win_mean = numpy.sum(numpy.sum(numpy.sum(img_strided, axis=-1), axis=-1), axis=-1) / (WIN_DIMS[0]*WIN_DIMS[1]*WIN_DIMS[2])
# As per http://en.wikipedia.org/wiki/Variance, we are removing the mean from every window,
# then squaring the result.
# Casting to 64 bit float inside, because the numbers (at least for our images) get pretty big
win_var = numpy.sum(numpy.sum(numpy.sum((( img_strided.T.astype('<f8') - win_mean.T.astype('<f8') )**2).T, axis=-1), axis=-1), axis=-1) / (WIN_DIMS[0]*WIN_DIMS[1]*WIN_DIMS[2])
# Prepare an output image of the right size, in order to replace the border removed with the windowed view call
out_img = numpy.zeros( DIMS, dtype='<f8' )
# copy borders out...
out_img[ WIN_DIMS[0]/2:DIMS[0]-WIN_DIMS[0]+1+WIN_DIMS[0]/2, WIN_DIMS[1]/2:DIMS[1]-WIN_DIMS[1]+1+WIN_DIMS[1]/2, WIN_DIMS[2]/2:DIMS[2]-WIN_DIMS[2]+1+WIN_DIMS[2]/2, ] = win_var
# output
return out_img.astype('>f4')
</code></pre>