Arbitrary n-dimensional rebin routine
![](https://secure.gravatar.com/avatar/ac3f41f2e11538be062ab6331fe68695.jpg?s=120&d=mm&r=g)
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@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.
![](https://secure.gravatar.com/avatar/ac3f41f2e11538be062ab6331fe68695.jpg?s=120&d=mm&r=g)
Okay, my first reply is to myself, with all the bits I forgot to add. Angus McMorland wrote:
[snip] import numpy as n import scipy.interpolate def nvals(n, dims): '''Returns xvals-like volumes, indexing over any dimension n. Will probably crash and burn if n >= len(dims).''' evList = '' for i in range( len(dims) - 1): if i < n: evList = evList + '%d' % dims[i] else: evList = evList + '%d' % dims[i + 1] if i < len(dims) - 2: evList = evList + ', ' evList = 'xvals( (%d,' % dims[n] + evList + ')' + ')' xs = eval( evList ) evList = '' for i in range( len(dims) ): if i < n: ind = i + 1 elif i == n: ind = 0 else: ind = i evList = evList + '%d' % ind if i < len(dims) - 1: evList = evList + ', ' evList = 'xs.transpose(' + evList + ')' return eval( evList )
def congrid(a, nud, method='neighbour', centre=False, minusone=False): '''Arbitrary resampling of source array to new dimension sizes.
[snip] Also the docstring should have something about the first couple of params: Usage: arg0: input array arg1: tuple of resulting dimensions Example: rebinned = rebin.congrid( raw, (2,2,2), \ method=[neighbour|linear|cubic|spline], \ minusone=[True|False], \ centre=[True|False]) Let's try again. -- Angus McMorland email a.mcmorland@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.
![](https://secure.gravatar.com/avatar/ac3f41f2e11538be062ab6331fe68695.jpg?s=120&d=mm&r=g)
Angus McMorland wrote:
Okay, my first reply is to myself, with all the bits I forgot to add.
Make that _some_ of the bits I forgot. Here's the last of them: def xvals(dims): '''Returns a volume with the appropriate x index value in each element.''' evList = ['n.fromfunction( xi, (%d' % dims[0] ] + \ [', %d' % dims[i+1] for i in range( len(dims) - 1 )] + \ [') )'] return eval( ''.join(evList) ) Apologies for the multiple postings. -- Angus McMorland email a.mcmorland@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.
![](https://secure.gravatar.com/avatar/764323a14e554c97ab74177e0bce51d4.jpg?s=120&d=mm&r=g)
Angus McMorland wrote:
Angus McMorland wrote:
Okay, my first reply is to myself, with all the bits I forgot to add.
Make that _some_ of the bits I forgot. Here's the last of them:
I'm not sure how many people are going to try to reconstruct the full module by copying and pasting from their email clients. Could you post it, in its entirety, to the SciPy wiki instead? This page would be good: http://www.scipy.org/Cookbook/Rebinning -- Robert Kern robert.kern@gmail.com "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
![](https://secure.gravatar.com/avatar/ac3f41f2e11538be062ab6331fe68695.jpg?s=120&d=mm&r=g)
Robert Kern wrote:
I'm not sure how many people are going to try to reconstruct the full module by copying and pasting from their email clients. Could you post it, in its entirety, to the SciPy wiki instead? This page would be good:
Done. Thanks for the advice, that's much better than email. The code is now up, in its entirety at the above page. A. -- Angus McMorland email a.mcmorland@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.
participants (2)
-
Angus McMorland
-
Robert Kern