search

PyWavelets / pywt

PyWavelets - Wavelet Transforms in Python
http://pywavelets.readthedocs.org
MIT License
2.08k stars 487 forks source link

Dual-tree complex wavelet transform #425

Open hgomersall opened 6 years ago

hgomersall commented 6 years ago

It just occurred to me that I wrote a pure python unencumbered implementation of the DT-CWT which might fit into this toolbox here.

It's written mostly for algorithmic understanding, not speed (so is very well documented/commented and cleanly implemented), but it is still usable. It was always my intention to write an optimised implementation at some point.

I can re-license if it makes sense.

grlee77 commented 6 years ago

Thanks @hgomersall. There was some prior discussion of adding this transform in #58 and #253.

Do you happen to recall the source of the q-shift filter coefficients in your implementation? I had previously used the ones from Dr. Selesnick's toolbox. It had filters of length 6, 14 and 18, but these were only specified to 8 decimal places (I think this was the precision in the table in one of the original publications).

Some time ago I had wanted to apply the DTCWT to some 4D data for a research project and was unable to find an implementation in dimensions > 3. I ended up developing an n-dimensional implementation, but haven't incorporated it into PyWavelets yet. It is in Python and was based on extending the 2D and 3D Matlab implementations by Dr. Selesnick that are available at http://eeweb.poly.edu/iselesni/WaveletSoftware/. Dr. Selesnick said he was fine with making a Python port based off of that code available under MIT license terms. I already had an option to use PyWavelets for the backend, I just haven't had a chance to clean it up and make a PR based on it yet.

There are also a couple of other Python-based implementations at: https://github.com/LaurentRDC/dualtree (1D, 2D) https://github.com/rjw57/dtcwt (1D, 2D, 3D)

hgomersall commented 6 years ago

The q-shifts I think I acquired from Nick Kingsbury, but they're published in one of his papers (though I'm not sure to quite the same precision).

It does sound like you might have it well covered. I'd be interested in your approach for 4D - the problem is finding a suitable data structure about 2D (frankly, even 2D isn't ideal!). I wasn't aware Ivan Selesnick had a toolbox, but that's good info.

Rich's (@rjw57) code is a direct translation of Nick Kingsbury's original matlab toolbox so cannot really be under a useful license. I wasn't aware of @LaurentRDC's code.

The original motivation for my code (which took me a while to publish) was to gain a deep understanding, so I'm sure it has use in the community, though I suspect not in PyWavelets. It also has an extension to complex input signals (which removes half space restriction on the outputs - real inputs generate outputs in which the negative frequency parts are the complex conjugates of the positive frequency parts, like the fourier transform). The complex input stuff was touched on in the IEEE DTCWT tutorial paper, but not really elaborated properly.

grlee77 commented 6 years ago

The original motivation for my code (which took me a while to publish) was to gain a deep understanding, so I'm sure it has use in the community, though I suspect not in PyWavelets.

I agree it definitely has use. Certainly, a 2D implementation is going to be easier to parse and understand than a generalized n-d implementation. Another bonus is that yours is already publicly available :)

I'd be interested in your approach for 4D - the problem is finding a suitable data structure about 2D (frankly, even 2D isn't ideal!)

The coefficients get stored in a manner analogous to wavedecn in PyWavelets. e.g. in PyWavelets, the set of detail coefficients corresponding to lowpass along axis 0 and highpass along axis 1 at a given level of decomposition would be indexed as coeffs[-level]['ad']. In the dual-tree case I use a tuple such as (0, 1) that is equal in length to the number of transformed axes and where 0 or 1 indicates which wavelet (e.g. 1 = dual wavelet) was used. So in 2D there are 4 sets of 'ad' coefficients at a given decomposition level indexed as:

coeffs[-level][(0, 0)]['ad'] coeffs[-level][(0, 1)]['ad'] coeffs[-level][(1, 0)]['ad'] coeffs[-level][(1, 1)]['ad'] In practice there is a linear recombination of the (0, 0), (1, 0), (0, 1), and (1, 1) trees to give the directional sensitivity as in the IEEE paper. I stored the recombined coefficients in the exact same structure, but this has the downside that now the (0, 0), (1, 1) are actually a linear combination and do not directly correspond to filtering with those filter combinations. It is also not obvious from that indexing scheme which element is sensitive to which direction.

This tuple-based key scheme above can be extended in a straightforward manner to n dimensions. The other detail is working out what the various weights are for the linear combinations to produce the desired real or complex directionally-selective transform. This has been published for 3D in another paper by Selesnick's group, but becomes a pain to work out by hand in high dimensions. I ended up precomputing for dimensions up to 6 using a SymPy script and storing those for later use to avoid a SymPy dependency at run time. It can in principle be computed up to any number of dimensions, but I have not used it on any real-world data in dimensions > 4.

tsikes commented 3 years ago

Has there been any progress on implementing this?