[SciPy-user] multidimensional wavelet packages

[Karl Young]

Hi Filip,

Thanks much (and thanks for the original package); I will go through the 
code and let you know if I come up with anything that would be worth 
incorporating (or let you know that your suggested addition works fine 
and should be added as is).

>Hi Karl,
>
>On Wed, Jan 7, 2009 at 22:15, Karl Young <Karl.Young at 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
>
><code>
>#!/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
>    characteristic.
>    """
>    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'))
>    else:
>        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')
>    pprint.pprint(result)
>
></code>
>
>Filip Wasilewski
>  
>


-- 

Karl Young
Center for Imaging of Neurodegenerative Diseases, UCSF          
VA Medical Center (114M)              Phone:  (415) 221-4810 x3114  lab        
4150 Clement Street                   FAX:    (415) 668-2864
San Francisco, CA 94121               Email:  karl young at ucsf edu