Aluode/PerceptionLabPortable
0
1import warnings
2from itertools import product
3
4import numpy as np
5
6from ._c99_config import _have_c99_complex
7from ._extensions._dwt import idwt_single
8from ._extensions._pywt import Modes, Wavelet, _check_dtype
9from ._extensions._swt import swt as _swt
10from ._extensions._swt import swt_axis as _swt_axis
11from ._extensions._swt import swt_max_level
12from ._multidim import idwt2, idwtn
13from ._utils import AxisError, _as_wavelet, _wavelets_per_axis
14
15__all__ = ["swt", "swt_max_level", 'iswt', 'swt2', 'iswt2', 'swtn', 'iswtn']
16
17
18def _rescale_wavelet_filterbank(wavelet, sf):
19 wav = Wavelet(wavelet.name + 'r',
20 [np.asarray(f) * sf for f in wavelet.filter_bank])
21
22 # copy attributes from the original wavelet
23 wav.orthogonal = wavelet.orthogonal
24 wav.biorthogonal = wavelet.biorthogonal
25 return wav
26
27
28def swt(data, wavelet, level=None, start_level=0, axis=-1,
29 trim_approx=False, norm=False):
30 """
31 Multilevel 1D stationary wavelet transform.
32
33 Parameters
34 ----------
35 data :
36 Input signal
37 wavelet :
38 Wavelet to use (Wavelet object or name)
39 level : int, optional
40 The number of decomposition steps to perform.
41 start_level : int, optional
42 The level at which the decomposition will begin (it allows one to
43 skip a given number of transform steps and compute
44 coefficients starting from start_level) (default: 0)
45 axis: int, optional
46 Axis over which to compute the SWT. If not given, the
47 last axis is used.
48 trim_approx : bool, optional
49 If True, approximation coefficients at the final level are retained.
50 norm : bool, optional
51 If True, transform is normalized so that the energy of the coefficients
52 will be equal to the energy of ``data``. In other words,
53 ``np.linalg.norm(data.ravel())`` will equal the norm of the
54 concatenated transform coefficients when ``trim_approx`` is True.
55
56 Returns
57 -------
58 coeffs : list
59 List of approximation and details coefficients pairs in order
60 similar to wavedec function::
61
62 [(cAn, cDn), ..., (cA2, cD2), (cA1, cD1)]
63
64 where n equals input parameter ``level``.
65
66 If ``start_level = m`` is given, then the beginning m steps are
67 skipped::
68
69 [(cAm+n, cDm+n), ..., (cAm+1, cDm+1), (cAm, cDm)]
70
71 If ``trim_approx`` is ``True``, then the output list is exactly as in
72 ``pywt.wavedec``, where the first coefficient in the list is the
73 approximation coefficient at the final level and the rest are the
74 detail coefficients::
75
76 [cAn, cDn, ..., cD2, cD1]
77
78 Notes
79 -----
80 The implementation here follows the "algorithm a-trous" and requires that
81 the signal length along the transformed axis be a multiple of ``2**level``.
82 If this is not the case, the user should pad up to an appropriate size
83 using a function such as ``numpy.pad``.
84
85 A primary benefit of this transform in comparison to its decimated
86 counterpart (``pywt.wavedecn``), is that it is shift-invariant. This comes
87 at cost of redundancy in the transform (the size of the output coefficients
88 is larger than the input).
89
90 When the following three conditions are true:
91
92 1. The wavelet is orthogonal
93 2. ``swt`` is called with ``norm=True``
94 3. ``swt`` is called with ``trim_approx=True``
95
96 the transform has the following additional properties that may be
97 desirable in applications:
98
99 1. energy is conserved
100 2. variance is partitioned across scales
101
102 When used with ``norm=True``, this transform is closely related to the
103 multiple-overlap DWT (MODWT) as popularized for time-series analysis,
104 although the underlying implementation is slightly different from the one
105 published in [1]_. Specifically, the implementation used here requires a
106 signal that is a multiple of ``2**level`` in length.
107
108 References
109 ----------
110 .. [1] DB Percival and AT Walden. Wavelet Methods for Time Series Analysis.
111 Cambridge University Press, 2000.
112 """
113
114 if not _have_c99_complex and np.iscomplexobj(data):
115 data = np.asarray(data)
116 kwargs = {"wavelet": wavelet, "level": level, "start_level": start_level,
117 "trim_approx": trim_approx, "axis": axis, "norm": norm}
118 coeffs_real = swt(data.real, **kwargs)
119 coeffs_imag = swt(data.imag, **kwargs)
120 if not trim_approx:
121 coeffs_cplx = []
122 for (cA_r, cD_r), (cA_i, cD_i) in zip(coeffs_real, coeffs_imag):
123 coeffs_cplx.append((cA_r + 1j*cA_i, cD_r + 1j*cD_i))
124 else:
125 coeffs_cplx = [cr + 1j*ci
126 for (cr, ci) in zip(coeffs_real, coeffs_imag)]
127 return coeffs_cplx
128
129 # accept array_like input; make a copy to ensure a contiguous array
130 dt = _check_dtype(data)
131 data = np.array(data, dtype=dt)
132
133 wavelet = _as_wavelet(wavelet)
134 if norm:
135 if not wavelet.orthogonal:
136 warnings.warn(
137 "norm=True, but the wavelet is not orthogonal: \n"
138 "\tThe conditions for energy preservation are not satisfied.")
139 wavelet = _rescale_wavelet_filterbank(wavelet, 1/np.sqrt(2))
140
141 if axis < 0:
142 axis = axis + data.ndim
143 if not 0 <= axis < data.ndim:
144 raise AxisError("Axis greater than data dimensions")
145
146 if level is None:
147 level = swt_max_level(data.shape[axis])
148
149 if data.ndim == 1:
150 ret = _swt(data, wavelet, level, start_level, trim_approx)
151 else:
152 ret = _swt_axis(data, wavelet, level, start_level, axis, trim_approx)
153 return ret
154
155
156def iswt(coeffs, wavelet, norm=False, axis=-1):
157 """
158 Multilevel 1D inverse discrete stationary wavelet transform.
159
160 Parameters
161 ----------
162 coeffs : array_like
163 Coefficients list of tuples::
164
165 [(cAn, cDn), ..., (cA2, cD2), (cA1, cD1)]
166
167 where cA is approximation, cD is details. Index 1 corresponds to
168 ``start_level`` from ``pywt.swt``.
169 wavelet : Wavelet object or name string
170 Wavelet to use
171 norm : bool, optional
172 Controls the normalization used by the inverse transform. This must
173 be set equal to the value that was used by ``pywt.swt`` to preserve the
174 energy of a round-trip transform.
175
176 Returns
177 -------
178 1D array of reconstructed data.
179
180 Examples
181 --------
182 >>> import pywt
183 >>> coeffs = pywt.swt([1,2,3,4,5,6,7,8], 'db2', level=2)
184 >>> pywt.iswt(coeffs, 'db2')
185 array([ 1., 2., 3., 4., 5., 6., 7., 8.])
186 """
187 # copy to avoid modification of input data
188 # If swt was called with trim_approx=False, first element is a tuple
189 trim_approx = not isinstance(coeffs[0], (tuple, list))
190 cA = coeffs[0] if trim_approx else coeffs[0][0]
191 if cA.ndim > 1:
192 # convert to swtn coefficient format and call iswtn
193 if trim_approx:
194 coeffs_nd = [cA] + [{'d': d} for d in coeffs[1:]]
195 else:
196 coeffs_nd = [{'a': a, 'd': d} for a, d in coeffs]
197 return iswtn(coeffs_nd, wavelet, axes=(axis,), norm=norm)
198 elif axis != 0 and axis != -1:
199 raise AxisError("Axis greater than data dimensions")
200 if not _have_c99_complex and np.iscomplexobj(cA):
201 if trim_approx:
202 coeffs_real = [c.real for c in coeffs]
203 coeffs_imag = [c.imag for c in coeffs]
204 else:
205 coeffs_real = [(ca.real, cd.real) for ca, cd in coeffs]
206 coeffs_imag = [(ca.imag, cd.imag) for ca, cd in coeffs]
207 kwargs = {"wavelet": wavelet, "norm": norm}
208 y = iswt(coeffs_real, **kwargs)
209 return y + 1j * iswt(coeffs_imag, **kwargs)
210
211 if trim_approx:
212 coeffs = coeffs[1:]
213
214 if cA.ndim != 1:
215 raise ValueError("iswt only supports 1D data")
216
217 dt = _check_dtype(cA)
218 output = np.array(cA, dtype=dt, copy=True)
219
220 # num_levels, equivalent to the decomposition level, n
221 num_levels = len(coeffs)
222 wavelet = _as_wavelet(wavelet)
223 if norm:
224 wavelet = _rescale_wavelet_filterbank(wavelet, np.sqrt(2))
225 mode = Modes.from_object('periodization')
226 for j in range(num_levels, 0, -1):
227 step_size = int(pow(2, j-1))
228 last_index = step_size
229 if trim_approx:
230 cD = coeffs[-j]
231 else:
232 _, cD = coeffs[-j]
233 cD = np.asarray(cD, dtype=_check_dtype(cD))
234 if cD.dtype != output.dtype:
235 # upcast to a common dtype (float64 or complex128)
236 if output.dtype.kind == 'c' or cD.dtype.kind == 'c':
237 dtype = np.complex128
238 else:
239 dtype = np.float64
240 output = np.asarray(output, dtype=dtype)
241 cD = np.asarray(cD, dtype=dtype)
242 for first in range(last_index): # 0 to last_index - 1
243
244 # Getting the indices that we will transform
245 indices = np.arange(first, len(cD), step_size)
246
247 # select the even indices
248 even_indices = indices[0::2]
249 # select the odd indices
250 odd_indices = indices[1::2]
251
252 # perform the inverse dwt on the selected indices,
253 # making sure to use periodic boundary conditions
254 # Note: indexing with an array of ints returns a contiguous
255 # copy as required by idwt_single.
256 x1 = idwt_single(output[even_indices],
257 cD[even_indices],
258 wavelet, mode)
259 x2 = idwt_single(output[odd_indices],
260 cD[odd_indices],
261 wavelet, mode)
262
263 # perform a circular shift right
264 x2 = np.roll(x2, 1)
265
266 # average and insert into the correct indices
267 output[indices] = (x1 + x2)/2.
268
269 return output
270
271
272def swt2(data, wavelet, level, start_level=0, axes=(-2, -1),
273 trim_approx=False, norm=False):
274 """
275 Multilevel 2D stationary wavelet transform.
276
277 Parameters
278 ----------
279 data : array_like
280 2D array with input data
281 wavelet : Wavelet object or name string, or 2-tuple of wavelets
282 Wavelet to use. This can also be a tuple of wavelets to apply per
283 axis in ``axes``.
284 level : int
285 The number of decomposition steps to perform.
286 start_level : int, optional
287 The level at which the decomposition will start (default: 0)
288 axes : 2-tuple of ints, optional
289 Axes over which to compute the SWT. Repeated elements are not allowed.
290 trim_approx : bool, optional
291 If True, approximation coefficients at the final level are retained.
292 norm : bool, optional
293 If True, transform is normalized so that the energy of the coefficients
294 will be equal to the energy of ``data``. In other words,
295 ``np.linalg.norm(data.ravel())`` will equal the norm of the
296 concatenated transform coefficients when ``trim_approx`` is True.
297
298 Returns
299 -------
300 coeffs : list
301 Approximation and details coefficients (for ``start_level = m``).
302 If ``trim_approx`` is ``False``, approximation coefficients are
303 retained for all levels::
304
305 [
306 (cA_m+level,
307 (cH_m+level, cV_m+level, cD_m+level)
308 ),
309 ...,
310 (cA_m+1,
311 (cH_m+1, cV_m+1, cD_m+1)
312 ),
313 (cA_m,
314 (cH_m, cV_m, cD_m)
315 )
316 ]
317
318 where cA is approximation, cH is horizontal details, cV is
319 vertical details, cD is diagonal details and m is ``start_level``.
320
321 If ``trim_approx`` is ``True``, approximation coefficients are only
322 retained at the final level of decomposition. This matches the format
323 used by ``pywt.wavedec2``::
324
325 [
326 cA_m+level,
327 (cH_m+level, cV_m+level, cD_m+level),
328 ...,
329 (cH_m+1, cV_m+1, cD_m+1),
330 (cH_m, cV_m, cD_m),
331 ]
332
333 Notes
334 -----
335 The implementation here follows the "algorithm a-trous" and requires that
336 the signal length along the transformed axes be a multiple of ``2**level``.
337 If this is not the case, the user should pad up to an appropriate size
338 using a function such as ``numpy.pad``.
339
340 A primary benefit of this transform in comparison to its decimated
341 counterpart (``pywt.wavedecn``), is that it is shift-invariant. This comes
342 at cost of redundancy in the transform (the size of the output coefficients
343 is larger than the input).
344
345 When the following three conditions are true:
346
347 1. The wavelet is orthogonal
348 2. ``swt2`` is called with ``norm=True``
349 3. ``swt2`` is called with ``trim_approx=True``
350
351 the transform has the following additional properties that may be
352 desirable in applications:
353
354 1. energy is conserved
355 2. variance is partitioned across scales
356
357 """
358 axes = tuple(axes)
359 data = np.asarray(data)
360 if len(axes) != 2:
361 raise ValueError("Expected 2 axes")
362 if len(axes) != len(set(axes)):
363 raise ValueError("The axes passed to swt2 must be unique.")
364 if data.ndim < len(np.unique(axes)):
365 raise ValueError("Input array has fewer dimensions than the specified "
366 "axes")
367
368 coefs = swtn(data, wavelet, level, start_level, axes, trim_approx, norm)
369 ret = []
370 if trim_approx:
371 ret.append(coefs[0])
372 coefs = coefs[1:]
373 for c in coefs:
374 if trim_approx:
375 ret.append((c['da'], c['ad'], c['dd']))
376 else:
377 ret.append((c['aa'], (c['da'], c['ad'], c['dd'])))
378 return ret
379
380
381def iswt2(coeffs, wavelet, norm=False, axes=(-2, -1)):
382 """
383 Multilevel 2D inverse discrete stationary wavelet transform.
384
385 Parameters
386 ----------
387 coeffs : list
388 Approximation and details coefficients::
389
390 [
391 (cA_n,
392 (cH_n, cV_n, cD_n)
393 ),
394 ...,
395 (cA_2,
396 (cH_2, cV_2, cD_2)
397 ),
398 (cA_1,
399 (cH_1, cV_1, cD_1)
400 )
401 ]
402
403 where cA is approximation, cH is horizontal details, cV is
404 vertical details, cD is diagonal details and n is the number of
405 levels. Index 1 corresponds to ``start_level`` from ``pywt.swt2``.
406 wavelet : Wavelet object or name string, or 2-tuple of wavelets
407 Wavelet to use. This can also be a 2-tuple of wavelets to apply per
408 axis.
409 norm : bool, optional
410 Controls the normalization used by the inverse transform. This must
411 be set equal to the value that was used by ``pywt.swt2`` to preserve
412 the energy of a round-trip transform.
413
414 Returns
415 -------
416 2D array of reconstructed data.
417
418 Examples
419 --------
420 >>> import pywt
421 >>> coeffs = pywt.swt2([[1,2,3,4],[5,6,7,8],
422 ... [9,10,11,12],[13,14,15,16]],
423 ... 'db1', level=2)
424 >>> pywt.iswt2(coeffs, 'db1')
425 array([[ 1., 2., 3., 4.],
426 [ 5., 6., 7., 8.],
427 [ 9., 10., 11., 12.],
428 [ 13., 14., 15., 16.]])
429
430 """
431
432 # If swt was called with trim_approx=False, first element is a tuple
433 trim_approx = not isinstance(coeffs[0], (tuple, list))
434 cA = coeffs[0] if trim_approx else coeffs[0][0]
435 if cA.ndim != 2 or axes != (-2, -1):
436 # convert to swtn coefficient format and call iswtn instead
437 if trim_approx:
438 coeffs_nd = [cA] + [{'da': h, 'ad': v, 'dd': d}
439 for h, v, d in coeffs[1:]]
440 else:
441 coeffs_nd = [{'aa': a, 'da': h, 'ad': v, 'dd': d}
442 for a, (h, v, d) in coeffs]
443 return iswtn(coeffs_nd, wavelet, axes=axes, norm=norm)
444 if not _have_c99_complex and np.iscomplexobj(cA):
445 if trim_approx:
446 coeffs_real = [cA.real]
447 coeffs_real += [(h.real, v.real, d.real) for h, v, d in coeffs[1:]]
448 coeffs_imag = [cA.imag]
449 coeffs_imag += [(h.imag, v.imag, d.imag) for h, v, d in coeffs[1:]]
450 else:
451 coeffs_real = [(a.real, (h.real, v.real, d.real))
452 for a, (h, v, d) in coeffs]
453 coeffs_imag = [(a.imag, (h.imag, v.imag, d.imag))
454 for a, (h, v, d) in coeffs]
455 kwargs = {"wavelet": wavelet, "norm": norm}
456 y = iswt2(coeffs_real, **kwargs)
457 return y + 1j * iswt2(coeffs_imag, **kwargs)
458
459 if trim_approx:
460 coeffs = coeffs[1:]
461
462 # copy to avoid modification of input data
463 dt = _check_dtype(cA)
464 output = np.array(cA, dtype=dt, copy=True)
465
466 if output.ndim != 2:
467 raise ValueError(
468 "iswt2 only supports 2D arrays. see iswtn for a general "
469 "n-dimensionsal ISWT")
470 # num_levels, equivalent to the decomposition level, n
471 num_levels = len(coeffs)
472 wavelets = _wavelets_per_axis(wavelet, axes=(0, 1))
473 if norm:
474 wavelets = [_rescale_wavelet_filterbank(wav, np.sqrt(2))
475 for wav in wavelets]
476
477 for j in range(num_levels):
478 step_size = int(pow(2, num_levels-j-1))
479 last_index = step_size
480 if trim_approx:
481 (cH, cV, cD) = coeffs[j]
482 else:
483 _, (cH, cV, cD) = coeffs[j]
484 # We are going to assume cH, cV, and cD are of equal size
485 if (cH.shape != cV.shape) or (cH.shape != cD.shape):
486 raise RuntimeError(
487 "Mismatch in shape of intermediate coefficient arrays")
488
489 # make sure output shares the common dtype
490 # (conversion of dtype for individual coeffs is handled within idwt2 )
491 common_dtype = np.result_type(*(
492 [dt, ] + [_check_dtype(c) for c in [cH, cV, cD]]))
493 if output.dtype != common_dtype:
494 output = output.astype(common_dtype)
495
496 for first_h in range(last_index): # 0 to last_index - 1
497 for first_w in range(last_index): # 0 to last_index - 1
498 # Getting the indices that we will transform
499 indices_h = slice(first_h, cH.shape[0], step_size)
500 indices_w = slice(first_w, cH.shape[1], step_size)
501
502 even_idx_h = slice(first_h, cH.shape[0], 2*step_size)
503 even_idx_w = slice(first_w, cH.shape[1], 2*step_size)
504 odd_idx_h = slice(first_h + step_size, cH.shape[0], 2*step_size)
505 odd_idx_w = slice(first_w + step_size, cH.shape[1], 2*step_size)
506
507 # perform the inverse dwt on the selected indices,
508 # making sure to use periodic boundary conditions
509 x1 = idwt2((output[even_idx_h, even_idx_w],
510 (cH[even_idx_h, even_idx_w],
511 cV[even_idx_h, even_idx_w],
512 cD[even_idx_h, even_idx_w])),
513 wavelets, 'periodization')
514 x2 = idwt2((output[even_idx_h, odd_idx_w],
515 (cH[even_idx_h, odd_idx_w],
516 cV[even_idx_h, odd_idx_w],
517 cD[even_idx_h, odd_idx_w])),
518 wavelets, 'periodization')
519 x3 = idwt2((output[odd_idx_h, even_idx_w],
520 (cH[odd_idx_h, even_idx_w],
521 cV[odd_idx_h, even_idx_w],
522 cD[odd_idx_h, even_idx_w])),
523 wavelets, 'periodization')
524 x4 = idwt2((output[odd_idx_h, odd_idx_w],
525 (cH[odd_idx_h, odd_idx_w],
526 cV[odd_idx_h, odd_idx_w],
527 cD[odd_idx_h, odd_idx_w])),
528 wavelets, 'periodization')
529
530 # perform a circular shifts
531 x2 = np.roll(x2, 1, axis=1)
532 x3 = np.roll(x3, 1, axis=0)
533 x4 = np.roll(x4, 1, axis=0)
534 x4 = np.roll(x4, 1, axis=1)
535 output[indices_h, indices_w] = (x1 + x2 + x3 + x4) / 4
536
537 return output
538
539
540def swtn(data, wavelet, level, start_level=0, axes=None, trim_approx=False,
541 norm=False):
542 """
543 n-dimensional stationary wavelet transform.
544
545 Parameters
546 ----------
547 data : array_like
548 n-dimensional array with input data.
549 wavelet : Wavelet object or name string, or tuple of wavelets
550 Wavelet to use. This can also be a tuple of wavelets to apply per
551 axis in ``axes``.
552 level : int
553 The number of decomposition steps to perform.
554 start_level : int, optional
555 The level at which the decomposition will start (default: 0)
556 axes : sequence of ints, optional
557 Axes over which to compute the SWT. A value of ``None`` (the
558 default) selects all axes. Axes may not be repeated.
559 trim_approx : bool, optional
560 If True, approximation coefficients at the final level are retained.
561 norm : bool, optional
562 If True, transform is normalized so that the energy of the coefficients
563 will be equal to the energy of ``data``. In other words,
564 ``np.linalg.norm(data.ravel())`` will equal the norm of the
565 concatenated transform coefficients when ``trim_approx`` is True.
566
567 Returns
568 -------
569 [{coeffs_level_n}, ..., {coeffs_level_1}]: list of dict
570 Results for each level are arranged in a dictionary, where the key
571 specifies the transform type on each dimension and value is a
572 n-dimensional coefficients array.
573
574 For example, for a 2D case the result at a given level will look
575 something like this::
576
577 {'aa': <coeffs> # A(LL) - approx. on 1st dim, approx. on 2nd dim
578 'ad': <coeffs> # V(LH) - approx. on 1st dim, det. on 2nd dim
579 'da': <coeffs> # H(HL) - det. on 1st dim, approx. on 2nd dim
580 'dd': <coeffs> # D(HH) - det. on 1st dim, det. on 2nd dim
581 }
582
583 For user-specified ``axes``, the order of the characters in the
584 dictionary keys map to the specified ``axes``.
585
586 If ``trim_approx`` is ``True``, the first element of the list contains
587 the array of approximation coefficients from the final level of
588 decomposition, while the remaining coefficient dictionaries contain
589 only detail coefficients. This matches the behavior of `pywt.wavedecn`.
590
591 Notes
592 -----
593 The implementation here follows the "algorithm a-trous" and requires that
594 the signal length along the transformed axes be a multiple of ``2**level``.
595 If this is not the case, the user should pad up to an appropriate size
596 using a function such as ``numpy.pad``.
597
598 A primary benefit of this transform in comparison to its decimated
599 counterpart (``pywt.wavedecn``), is that it is shift-invariant. This comes
600 at cost of redundancy in the transform (the size of the output coefficients
601 is larger than the input).
602
603 When the following three conditions are true:
604
605 1. The wavelet is orthogonal
606 2. ``swtn`` is called with ``norm=True``
607 3. ``swtn`` is called with ``trim_approx=True``
608
609 the transform has the following additional properties that may be
610 desirable in applications:
611
612 1. energy is conserved
613 2. variance is partitioned across scales
614
615 """
616 data = np.asarray(data)
617 if not _have_c99_complex and np.iscomplexobj(data):
618 kwargs = {"wavelet": wavelet, "level": level, "start_level": start_level,
619 "trim_approx": trim_approx, "axes": axes, "norm": norm}
620 real = swtn(data.real, **kwargs)
621 imag = swtn(data.imag, **kwargs)
622 if trim_approx:
623 cplx = [real[0] + 1j * imag[0]]
624 offset = 1
625 else:
626 cplx = []
627 offset = 0
628 for rdict, idict in zip(real[offset:], imag[offset:]):
629 cplx.append(
630 {k: rdict[k] + 1j * idict[k] for k in rdict})
631 return cplx
632
633 if data.dtype == np.dtype('object'):
634 raise TypeError("Input must be a numeric array-like")
635 if data.ndim < 1:
636 raise ValueError("Input data must be at least 1D")
637
638 if axes is None:
639 axes = range(data.ndim)
640 axes = [a + data.ndim if a < 0 else a for a in axes]
641 if any(a < 0 or a >= data.ndim for a in axes):
642 raise AxisError("Axis greater than data dimensions")
643 if len(axes) != len(set(axes)):
644 raise ValueError("The axes passed to swtn must be unique.")
645 num_axes = len(axes)
646
647 wavelets = _wavelets_per_axis(wavelet, axes)
648 if norm:
649 if not np.all([wav.orthogonal for wav in wavelets]):
650 warnings.warn(
651 "norm=True, but the wavelets used are not orthogonal: \n"
652 "\tThe conditions for energy preservation are not satisfied.")
653 wavelets = [_rescale_wavelet_filterbank(wav, 1/np.sqrt(2))
654 for wav in wavelets]
655 ret = []
656 for i in range(start_level, start_level + level):
657 coeffs = [('', data)]
658 for axis, wavelet in zip(axes, wavelets):
659 new_coeffs = []
660 for subband, x in coeffs:
661 cA, cD = _swt_axis(x, wavelet, level=1, start_level=i,
662 axis=axis)[0]
663 new_coeffs.extend([(subband + 'a', cA),
664 (subband + 'd', cD)])
665 coeffs = new_coeffs
666
667 coeffs = dict(coeffs)
668 ret.append(coeffs)
669
670 # data for the next level is the approximation coeffs from this level
671 data = coeffs['a' * num_axes]
672 if trim_approx:
673 coeffs.pop('a' * num_axes)
674 if trim_approx:
675 ret.append(data)
676 ret.reverse()
677 return ret
678
679
680def iswtn(coeffs, wavelet, axes=None, norm=False):
681 """
682 Multilevel nD inverse discrete stationary wavelet transform.
683
684 Parameters
685 ----------
686 coeffs : list
687 [{coeffs_level_n}, ..., {coeffs_level_1}]: list of dict
688 wavelet : Wavelet object or name string, or tuple of wavelets
689 Wavelet to use. This can also be a tuple of wavelets to apply per
690 axis in ``axes``.
691 axes : sequence of ints, optional
692 Axes over which to compute the inverse SWT. Axes may not be repeated.
693 The default is ``None``, which means transform all axes
694 (``axes = range(data.ndim)``).
695 norm : bool, optional
696 Controls the normalization used by the inverse transform. This must
697 be set equal to the value that was used by ``pywt.swtn`` to preserve
698 the energy of a round-trip transform.
699
700 Returns
701 -------
702 nD array of reconstructed data.
703
704 Examples
705 --------
706 >>> import pywt
707 >>> coeffs = pywt.swtn([[1,2,3,4],[5,6,7,8],
708 ... [9,10,11,12],[13,14,15,16]],
709 ... 'db1', level=2)
710 >>> pywt.iswtn(coeffs, 'db1')
711 array([[ 1., 2., 3., 4.],
712 [ 5., 6., 7., 8.],
713 [ 9., 10., 11., 12.],
714 [ 13., 14., 15., 16.]])
715
716 """
717
718 # key length matches the number of axes transformed
719 ndim_transform = max(len(key) for key in coeffs[-1])
720 trim_approx = not isinstance(coeffs[0], dict)
721 cA = coeffs[0] if trim_approx else coeffs[0]['a'*ndim_transform]
722
723 if not _have_c99_complex and np.iscomplexobj(cA):
724 if trim_approx:
725 coeffs_real = [coeffs[0].real]
726 coeffs_imag = [coeffs[0].imag]
727 coeffs = coeffs[1:]
728 else:
729 coeffs_real = []
730 coeffs_imag = []
731 coeffs_real += [{k: v.real for k, v in c.items()} for c in coeffs]
732 coeffs_imag += [{k: v.imag for k, v in c.items()} for c in coeffs]
733 kwargs = {"wavelet": wavelet, "axes": axes, "norm": norm}
734 y = iswtn(coeffs_real, **kwargs)
735 return y + 1j * iswtn(coeffs_imag, **kwargs)
736
737 if trim_approx:
738 coeffs = coeffs[1:]
739
740 # copy to avoid modification of input data
741 dt = _check_dtype(cA)
742 output = np.array(cA, dtype=dt, copy=True)
743 ndim = output.ndim
744
745 if axes is None:
746 axes = range(output.ndim)
747 axes = [a + ndim if a < 0 else a for a in axes]
748 if len(axes) != len(set(axes)):
749 raise ValueError("The axes passed to swtn must be unique.")
750 if ndim_transform != len(axes):
751 raise ValueError("The number of axes used in iswtn must match the "
752 "number of dimensions transformed in swtn.")
753
754 # num_levels, equivalent to the decomposition level, n
755 num_levels = len(coeffs)
756 wavelets = _wavelets_per_axis(wavelet, axes)
757 if norm:
758 wavelets = [_rescale_wavelet_filterbank(wav, np.sqrt(2))
759 for wav in wavelets]
760
761 # initialize various slice objects used in the loops below
762 # these will remain slice(None) only on axes that aren't transformed
763 indices = [slice(None), ]*ndim
764 even_indices = [slice(None), ]*ndim
765 odd_indices = [slice(None), ]*ndim
766 odd_even_slices = [slice(None), ]*ndim
767
768 for j in range(num_levels):
769 step_size = int(pow(2, num_levels-j-1))
770 last_index = step_size
771 if not trim_approx:
772 a = coeffs[j].pop('a'*ndim_transform) # will restore later
773 details = coeffs[j]
774 # make sure dtype matches the coarsest level approximation coefficients
775 common_dtype = np.result_type(*(
776 [dt, ] + [v.dtype for v in details.values()]))
777 if output.dtype != common_dtype:
778 output = output.astype(common_dtype)
779
780 # We assume all coefficient arrays are of equal size
781 shapes = [v.shape for k, v in details.items()]
782 if len(set(shapes)) != 1:
783 raise RuntimeError(
784 "Mismatch in shape of intermediate coefficient arrays")
785
786 # shape of a single coefficient array, excluding non-transformed axes
787 coeff_trans_shape = tuple([shapes[0][ax] for ax in axes])
788
789 # nested loop over all combinations of axis offsets at this level
790 for firsts in product(*([range(last_index), ]*ndim_transform)):
791 for first, sh, ax in zip(firsts, coeff_trans_shape, axes):
792 indices[ax] = slice(first, sh, step_size)
793 even_indices[ax] = slice(first, sh, 2*step_size)
794 odd_indices[ax] = slice(first+step_size, sh, 2*step_size)
795
796 # nested loop over all combinations of odd/even inidices
797 approx = output.copy()
798 output[tuple(indices)] = 0
799 ntransforms = 0
800 for odds in product(*([(0, 1), ]*ndim_transform)):
801 for o, ax in zip(odds, axes):
802 if o:
803 odd_even_slices[ax] = odd_indices[ax]
804 else:
805 odd_even_slices[ax] = even_indices[ax]
806 # extract the odd/even indices for all detail coefficients
807 details_slice = {}
808 for key, value in details.items():
809 details_slice[key] = value[tuple(odd_even_slices)]
810 details_slice['a'*ndim_transform] = approx[
811 tuple(odd_even_slices)]
812
813 # perform the inverse dwt on the selected indices,
814 # making sure to use periodic boundary conditions
815 x = idwtn(details_slice, wavelets, 'periodization', axes=axes)
816 for o, ax in zip(odds, axes):
817 # circular shift along any odd indexed axis
818 if o:
819 x = np.roll(x, 1, axis=ax)
820 output[tuple(indices)] += x
821 ntransforms += 1
822 output[tuple(indices)] /= ntransforms # normalize
823 if not trim_approx:
824 coeffs[j]['a'*ndim_transform] = a # restore approx coeffs to dict
825 return output
826