Python pywt.waverec2() Examples

def generate_basis():
    """generate the basis"""
    x = np.zeros((56, 56))
    coefs = pywt.wavedec2(x, 'db1')
    n_levels = len(coefs)
    basis = []
    for i in range(n_levels):
        coefs[i] = list(coefs[i])
        n_filters = len(coefs[i])
        for j in range(n_filters):
            for m in range(coefs[i][j].shape[0]):
                try:
                    for n in range(coefs[i][j].shape[1]):
                        coefs[i][j][m][n] = 1
                        temp_basis = pywt.waverec2(coefs, 'db1')
                        basis.append(temp_basis)
                        coefs[i][j][m][n] = 0
                except IndexError:
                    coefs[i][j][m] = 1
                    temp_basis = pywt.waverec2(coefs, 'db1')
                    basis.append(temp_basis)
                    coefs[i][j][m] = 0

    basis = np.array(basis)
    return basis 

def hfilter(diff_image, var_image, threshold=1, ndamp=10):
    """
    This code was inspired from: https://github.com/spacetelescope/sprint_notebooks/blob/master/lucy_damped_haar.ipynb
    I believe it was initially written by Justin Ely: https://github.com/justincely
    It was buggy and not working properly with every image sizes.
    I have thus exchanged it by using pyWavelet (pywt) and a custom function htrans
    to calculate the matrix for the var_image.
    """
    him, coeff_slices = pywt.coeffs_to_array(pywt.wavedec2(diff_image.astype(np.float), 'haar'), padding=0)
    dvarim = htrans(var_image.astype(np.float))
    
    sqhim = ((him/threshold)**2)/dvarim
    index = np.where(sqhim < 1)
    
    if len(index[0]) == 0:
        return diff_image
    
    # Eq. 8 of White is derived leading to N*x^(N-1)-(N-1)*x^N  :DOI: 10.1117/12.176819
    sqhim = sqhim[index] * (ndamp * sqhim[index]**(ndamp-1) - (ndamp-1)*sqhim[index]**ndamp)
    him[index] = sign(threshold*np.sqrt(dvarim[index] * sqhim), him[index])
    
    return pywt.waverec2(pywt.array_to_coeffs(him, coeff_slices, output_format='wavedec2'), 'haar')[:diff_image.shape[0],:diff_image.shape[1]] 

def fusion(self):
        self._load_images()
        coeffss = []
        for image in self._images:
            coeffss.append(pywt.wavedec2(image, 'db1', level=self._zt))
        # low pass
        if self._mp == 0:
            cAF = coeffss[0][0]
            for coeffs in coeffss[1:]:
                cAF += coeffs[0]
            cAF = cAF/len(coeffs)
        # high pass
        if self._ap == 2:
            hipassF  = coeffss[0][1:]
            for coeffs in coeffss[1:]:   # every image
                for idxLevel, HVDs in enumerate(coeffs[1:]):   # every level
                    for idxDirec, HVD in enumerate(HVDs):
                        maxMap = hipassF[idxLevel][idxDirec] < HVD
                        hipassF[idxLevel][idxDirec][maxMap] = HVD[maxMap]

        coeffsFusion = [cAF,] + hipassF
        self._fusionImage = pywt.waverec2(coeffsFusion, 'db1')
        return self._fusionImage 

def test_wavedecn_coeff_reshape_even():
    # verify round trip is correct:
    #   wavedecn - >coeffs_to_array-> array_to_coeffs -> waverecn
    # This is done for wavedec{1, 2, n}
    rng = np.random.RandomState(1234)
    params = {'wavedec': {'d': 1, 'dec': pywt.wavedec, 'rec': pywt.waverec},
              'wavedec2': {'d': 2, 'dec': pywt.wavedec2, 'rec': pywt.waverec2},
              'wavedecn': {'d': 3, 'dec': pywt.wavedecn, 'rec': pywt.waverecn}}
    N = 28
    for f in params:
        x1 = rng.randn(*([N] * params[f]['d']))
        for mode in pywt.Modes.modes:
            for wave in wavelist:
                w = pywt.Wavelet(wave)
                maxlevel = pywt.dwt_max_level(np.min(x1.shape), w.dec_len)
                if maxlevel == 0:
                    continue

                coeffs = params[f]['dec'](x1, w, mode=mode)
                coeff_arr, coeff_slices = pywt.coeffs_to_array(coeffs)
                coeffs2 = pywt.array_to_coeffs(coeff_arr, coeff_slices,
                                               output_format=f)
                x1r = params[f]['rec'](coeffs2, w, mode=mode)

                assert_allclose(x1, x1r, rtol=1e-4, atol=1e-4) 

def test_waverec2_odd_length():
    x = np.ones((10, 6))
    coeffs = pywt.wavedec2(x, 'db1')
    assert_allclose(pywt.waverec2(coeffs, 'db1'), x, rtol=1e-12) 

def _wavelet_decomp(X, w):
    c = pywt.wavedec2(X, w)

    X_out = []
    for k in range(len(c)):
        c_ = c[:]
        for k_ in set(range(len(c))).difference([k]):
            c_[k_] = tuple([np.zeros_like(v) for v in c[k_]])
        X_k = pywt.waverec2(c_, w)
        X_out.append(X_k)

    return X_out 

def analyze(self, images):
        image = images['gray']

        if self.image_scale is not None:
            assert self.image_scale & (self.image_scale - 1) == 0, "image_scale is not power of 2"
        else:
            image_natural_scale = 2 ** int(numpy.log2(min(image.shape[0:2])))
            image_scale = max(image_natural_scale, self.hash_size)

        ll_max_level = int(numpy.log2(image_scale))

        level = int(numpy.log2(self.hash_size))
        assert self.hash_size & (self.hash_size - 1) == 0, "hash_size is not power of 2"
        assert level <= ll_max_level, "hash_size in a wrong range"
        dwt_level = ll_max_level - level

        scaled = cv2.resize(image, (image_scale, image_scale), interpolation=cv2.INTER_AREA)
        pixels = numpy.float32(scaled)
        pixels /= 255

        # Remove low level frequency LL(max_ll) if @remove_max_haar_ll using haar filter
        if self.remove_max_haar_ll:
            coeffs = pywt.wavedec2(pixels, 'haar', level=ll_max_level)
            coeffs = list(coeffs)
            coeffs[0] *= 0
            pixels = pywt.waverec2(coeffs, 'haar')

        # Use LL(K) as freq, where K is log2(@hash_size)
        coeffs = pywt.wavedec2(pixels, self.mode, level=dwt_level)
        dwt_low = coeffs[0]

        # Substract median and compute hash
        med = numpy.median(dwt_low)
        diff = dwt_low > med
        return numpy.packbits(numpy.uint8(diff.reshape(-1, 1))) 

def test_waverec2_axes_errors():
    data = np.ones((4, 4))
    c = pywt.wavedec2(data, 'haar')
    # integer axes not allowed
    assert_raises(TypeError, pywt.waverec2, c, 'haar', axes=1)
    # non-unique axes not allowed
    assert_raises(ValueError, pywt.waverec2, c, 'haar', axes=(0, 0))
    # out of range axis not allowed
    assert_raises(ValueError, pywt.waverec2, c, 'haar', axes=(0, 2)) 

def test_waverec2_axes_subsets():
    rstate = np.random.RandomState(0)
    data = rstate.standard_normal((8, 8, 8))
    # test all combinations of 2 out of 3 axes transformed
    for axes in combinations((0, 1, 2), 2):
        coefs = pywt.wavedec2(data, 'haar', axes=axes)
        rec = pywt.waverec2(coefs, 'haar', axes=axes)
        assert_allclose(rec, data, atol=1e-14) 

def test_waverec2_none_coeffs():
    x = np.arange(24).reshape(6, 4)
    coeffs = pywt.wavedec2(x, 'db1')
    coeffs[1] = (None, None, None)
    assert_(x.shape == pywt.waverec2(coeffs, 'db1').shape)

####
# nd multilevel dwt function tests
#### 

def recon(self,z1):
        """
        Wavelet reconstruction:  coefficients -> image
        """
        coeffs = pywt.array_to_coeffs(z1, self.coeff_slices, \
            output_format='wavedec2')
        z0 = pywt.waverec2(coeffs, wavelet=self.wavelet, mode=self.mode)
        return z0 

def test_waverec2_invalid_inputs():
    # input must be list or tuple
    assert_raises(ValueError, pywt.waverec2, np.ones((8, 8)), 'haar')

    # input list cannot be empty
    assert_raises(ValueError, pywt.waverec2, [], 'haar') 

def test_wavedec2_complex():
    data = np.ones((4, 4)) + 1j
    coeffs = pywt.wavedec2(data, 'db1')
    assert_(len(coeffs) == 3)
    assert_allclose(pywt.waverec2(coeffs, 'db1'), data, rtol=1e-12) 

def test_waverec2_all_wavelets_modes():
    # test 2D case using all wavelets and modes
    rstate = np.random.RandomState(1234)
    r = rstate.randn(80, 96)
    for wavelet in wavelist:
        for mode in pywt.Modes.modes:
            coeffs = pywt.wavedec2(r, wavelet, mode=mode)
            assert_allclose(pywt.waverec2(coeffs, wavelet, mode=mode),
                            r, rtol=tol_single, atol=tol_single) 

def test_multilevel_dtypes_2d():
    wavelet = pywt.Wavelet('haar')
    for dt_in, dt_out in zip(dtypes_in, dtypes_out):
        # wavedec2, waverec2
        x = np.ones((8, 8), dtype=dt_in)
        errmsg = "wrong dtype returned for {0} input".format(dt_in)
        cA, coeffsD2, coeffsD1 = pywt.wavedec2(x, wavelet, level=2)
        assert_(cA.dtype == dt_out, "wavedec2: " + errmsg)
        for c in coeffsD1:
            assert_(c.dtype == dt_out, "wavedec2: " + errmsg)
        for c in coeffsD2:
            assert_(c.dtype == dt_out, "wavedec2: " + errmsg)
        x_roundtrip = pywt.waverec2([cA, coeffsD2, coeffsD1], wavelet)
        assert_(x_roundtrip.dtype == dt_out, "waverec2: " + errmsg) 

def test_waverec2_accuracies():
    rstate = np.random.RandomState(1234)
    x0 = rstate.randn(4, 4)
    for dt, tol in dtypes_and_tolerances:
        x = x0.astype(dt)
        if np.iscomplexobj(x):
            x += 1j*rstate.randn(4, 4).astype(x.real.dtype)
        coeffs = pywt.wavedec2(x, 'db1')
        assert_(len(coeffs) == 3)
        assert_allclose(pywt.waverec2(coeffs, 'db1'), x, atol=tol, rtol=tol) 

def generate_basis():
    """generate the basis"""
    x = np.zeros((64, 64))
    coefs = pywt.wavedec2(x, 'db1')
    n_levels = len(coefs)
    basis = []
    for i in range(n_levels):
        coefs[i] = list(coefs[i])
        n_filters = len(coefs[i])
        for j in range(n_filters):
            for m in range(coefs[i][j].shape[0]):
                try:
                    for n in range(coefs[i][j].shape[1]):
                        coefs[i][j][m][n] = 1
                        temp_basis = pywt.waverec2(coefs, 'db1')
                        basis.append(temp_basis)
                        coefs[i][j][m][n] = 0
                except IndexError:
                    coefs[i][j][m] = 1
                    temp_basis = pywt.waverec2(coefs, 'db1')
                    basis.append(temp_basis)
                    coefs[i][j][m] = 0

    basis = np.array(basis)
    return basis 

def get_image(coefs_list):
    x = np.zeros((64, 64, 3))
    for i in range(3):
        x[:, :, i] = pywt.waverec2(coefs_list[i], 'db1')
    return x 

def _rmatvec(self, x):
        x = np.reshape(x, self.dimsd)
        x = pywt.array_to_coeffs(x, self.sl, output_format='wavedec2')
        y = pywt.waverec2(x, wavelet=self.waveletadj, mode='periodization',
                          axes=self.dirs)
        y = self.pad.rmatvec(y.ravel())
        return y