[SciPy-dev] Arbitrary n-dimensional rebin routine

[Angus McMorland]

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.