[SciPy-dev] Arbitrary n-dimensional rebin routine
Angus McMorland
a.mcmorland at auckland.ac.nz
Tue Apr 11 22:10:35 CDT 2006
Hi all,
I'm not sure of the protocol here, so I'll try the list first...
Here's an n-dimensional rebinning routine I've written (based loosely on
IDLs congrid function), and hope it might be of some use. Perhaps it can
get bundled with Gary Ruben's rebin function, currently on the wiki
(http://scipy.org/Cookbook/Rebinning). I've thrown a few simple tests at
the routine so far, with good results, but it could probably do with
more field testing.
The routine is currently optimised to take as long and as much memory as
possible, so if someone with more knowledge of these things than me is
interested to refine it a bit, that'd be great. As a scipy beginner, I'm
also really interested to get some feedback on coding best-practices
that I should have used and haven't.
Cheers,
Angus.
def congrid(a, nud, method='neighbour', centre=False, minusone=False):
'''Arbitrary resampling of source array to new dimension sizes.
Currently only supports maintaining the same number of dimensions.
Also doesn''t work for 1-D arrays unless promoted to shape (x,1).
Based loosely on IDL''s congrid routine, which apparently originally
came from a VAX/VMS routine of the same name.
I''m not completely sure of the validity of using parallel 1-D
interpolations repeated along each axis in succession, but the
results are visually compelling.
method:
neighbour - closest value from original data
the ''kinds'' supported by scipy.interpolate.interp1d
centre:
True - interpolation points are at the centres of the bins
False - points are at the front edge of the bin
minusone:
For example- inarray.shape = (i,j) & new dimensions = (x,y)
False - inarray is resampled by factors of (i/x) * (j/y)
True - inarray is resampled by(i-1)/(x-1) * (j-1)/(y-1)
This prevents extrapolation one element beyond bounds of input
array.
'''
if not a.dtype.type in [n.typeDict['Float32'], n.typeDict['Float64']]:
print "Converting to float"
a = a.astype('Float32')
if minusone:
m1 = 1.
else:
m1 = 0.
if centre:
ofs = 0.5
else:
ofs = 0.
old = n.asarray( a.shape )
ndims = len( old )
if len( nud ) != ndims:
print "Congrid: dimensions error. This routine currently only
support rebinning to the same number of dimensions."
return None
nudr = n.asarray( nud ).astype('Float32')
dimlist = []
if method == 'neighbour':
for i in range( ndims ):
base = nvals(i, nudr)
dimlist.append( (old[i] - m1) / (nudr[i] - m1) \
* (base + ofs) - ofs )
cd = n.array( dimlist )
cdr = cd.round().astype( 'UInt16' )
nua = a[list( cdr )]
return nua
elif method in ['nearest','linear','cubic','spline']:
# calculate new dims
for i in range( ndims ):
base = n.arange( nudr[i] )
dimlist.append( (old[i] - m1) / (nudr[i] - m1) \
* (base + ofs) - ofs )
# specify old dims
olddims = [n.arange(i).astype('Float32') for i in list( a.shape )]
# first interpolation - for ndims = any
mint = scipy.interpolate.interp1d( olddims[-1], a, kind=method )
nua = mint( dimlist[-1] )
trorder = [ndims - 1] + range( ndims - 1 )
for i in range( ndims - 2, -1, -1 ):
nua = nua.transpose( trorder )
mint = scipy.interpolate.interp1d( olddims[i], nua,
kind=method )
nua = mint( dimlist[i] )
if ndims > 1:
# need one more transpose to return to original dimensions
nua = nua.transpose( trorder )
return nua
else:
print "Congrid error: Unrecognized interpolation type.\n", \
"This routine currently only supports
\'nearest\',\'linear\',\'cubic\', and \'spline\'."
return None
--
Angus McMorland
email a.mcmorland at auckland.ac.nz
mobile +64-21-155-4906
PhD Student, Neurophysiology / Multiphoton & Confocal Imaging
Physiology, University of Auckland
phone +64-9-3737-599 x89707
Armourer, Auckland University Fencing
Secretary, Fencing North Inc.
More information about the Scipy-dev
mailing list