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.