[SciPy-user] multidimensional wavelet packages

Filip Wasilewski filipwasilewski@gmail....
Wed Jan 7 18:14:00 CST 2009

Hi Karl,

On Wed, Jan 7, 2009 at 22:15, Karl Young <Karl.Young@ucsf.edu> wrote:
> Hi Stefan,
> Thanks; I'd looked a little at PyWavelets and figured that what you
> suggest might be what I ended up hacking but thought maybe some
> enterprising neuroimager (or other person working with 3D, 4D data)
> might have already done so :-)

I haven't seen a 3D transform implementation in Python, but I can give
you some hints on extending PyWavelets.

First of all take a look at 2D DWT and IDWT implementation at [1]. It
follows a standard pattern of transforming rows and then columns using
1D transform and producing 2D arrays of coefficients.

As you may already know, to perform n-dimensional transform one will
have to apply 1D transform over each dimension, doubling the number of
n-dimensional coefficient arrays with every step (approximation and
details coefficients) -- see [2] for a 3D example.

Below is a naive implementation of this algorithm. As you can see it
is very short and even seems to work (I have verified results for 2d
case only), but unfortunately it has several major drawbacks in the
recursive approach (worst possible memory management and twice the
necessary computations because of PyWavelets missing true downcoef_a
and downcoef_d functions for use with apply_along_axis[3]).

I think it could be converted into something like the dwt2 from [1]
with freeing intermediate arrays, but I guess the resulting code may
become very complex, so the solution with apply_along_axis is still
very attractive (it only needs adding optimized downcoef_a and
downcoef_d functions to PyWavelets and converting recursion into
iteration to better handle memory usage).

Let me know if you come out with a more optimal solution, so if you
agree I could include it in PyWavelets.

Hope that will help you with n-dimensional implementation.

[1] http://projects.scipy.org/wavelets/browser/pywt/trunk/pywt/multidim.py
[2] http://taco.poly.edu/WaveletSoftware/standard3D.html
[3] http://docs.scipy.org/doc/numpy/reference/generated/numpy.apply_along_axis.html

#!/usr/bin/env python
# Author: Filip Wasilewski
# Licence: Public Domain

import numpy
import pywt

# Helpers for numpy.apply_along_axis, which expects a 1D array
# as the function output
def downcoef_a(*args, **kwargs):
    """Returns DWT approximation coeffs."""
    return pywt.dwt(*args, **kwargs)[0]

def downcoef_d(*args, **kwargs):
    """Returns DWT details coeffs."""
    return pywt.dwt(*args, **kwargs)[1]

def dwt_n(data, wavelet, mode='sym', axis=0, subband=''):
    """N-dimensional Discrete Wavelet Transform

    Note: This is a proof of concept with worst possible memory usage
    dim = len(data.shape)
    if axis < dim:
        cA = numpy.apply_along_axis(downcoef_a, axis, data, wavelet, mode)
        cD = numpy.apply_along_axis(downcoef_d, axis, data, wavelet, mode)
        return (dwt_n(cA, wavelet, mode, axis+1, subband=subband+'L'),
                dwt_n(cD, wavelet, mode, axis+1, subband=subband+'H'))
        return (subband, data) # (subband name, coeffs)

if __name__ == '__main__':
    import pprint
    x = numpy.ones((4, 4, 4, 4)) # 4D array
    result = dwt_n(x, 'db1')


Filip Wasilewski

>>Hi Karl
>>The only Python wavelet library I've ever used is the one at
>>It works pretty well, if only for 1D and 2D cases.  IIRC, some of the
>>orthogonal wavelet transforms are separable, so you may be able to
>>construct a 3D transform using the 1D functions already implemented.
>>2009/1/7 Young, Karl <karl.young@ucsf.edu>:
>>>I was just curious re. whether anyone on the list is aware of any multidimensional wavelet packages (either in python or with a python interface) - 3D and 4D is mainly what I'm looking for. I've searched a little and know there has been discussion and some development of wavelet packages for SciPy but I haven't kept up with that and it didn't look like their was anything multidimensional currently available or in the works. I don't need anything terribly fancy (e.g. just Haar wavelets would suffice) and can probably hack something but thought I'd check, both for selfish reasons (lazy !) and because if I'm going to do any work, contributing to a community effort should any already exist, is certainly preferable.

More information about the SciPy-user mailing list