def denoise():
    wave = 'db4'
    sig = 20
    tau1 = 3*sig
    tau2 = 3*sig/2
    noisyLena = lena + np.random.normal(scale = sig, size=lena.shape)
    lw = pywt.wavedec2(noisyLena, wave, level=4)
    lwt1 = hardThresh(lw, tau1)
    lwt2 = softThresh(lw, tau2)
    rlena1 = pywt.waverec2(lwt1, wave)
    rlena2 = pywt.waverec2(lwt2, wave)
    plt.subplot(131)
    plt.imshow(noisyLena, cmap=plt.cm.Greys_r)
    plt.axis('off')
    
    plt.subplot(132)
    plt.imshow(rlena1, cmap=plt.cm.Greys_r)
    plt.axis('off')
    
    plt.subplot(133)
    plt.imshow(rlena2, cmap=plt.cm.Greys_r)
    plt.axis('off')
    
    plt.savefig('denoise.pdf')
    plt.clf()

def create_haar_dictionary(p=8):
    import pywt
    c = pywt.wavedec2(np.zeros((p, p)), 'haar')
    D = []
    for k in range(1, len(c)):
        for i in range(3):
            ck = c[k][i]
            l = ck.shape[0]
            for j in range(l):
                for m in range(l):
                    ck[j, m] = 1
                    D += [pywt.waverec2(c, 'haar')]
                    ck[j, m] = 0
    ck = c[0]
    l = ck.shape[0]
    for j in range(l):
        for m in range(l):
            ck[j, m] = 1
            D += [pywt.waverec2(c, 'haar')]
            ck[j, m] = 0
    D = np.array(D).reshape(-1, p*p)

    Dn = []
    for i in range(15):
        Dn += _translate(D[i].reshape((p, p)))
    Dn = np.array(Dn).reshape((-1, p*p))
    i0 = np.sum(abs(Dn), axis=1) != 0
    return Dn[i0]

def w2dReg(img, mode=mode, level=level, noiseSigma = noiseSigma):
    # compute coefficients
    noisy_img = add_noise(img, noiseSigma = noiseSigma)
    rec_coeffs = denoise(noisy_img, mode, level, noiseSigma)
    # reconstruction
    rec_img = pywt.waverec2(rec_coeffs, mode);
    return noisy_img, rec_img, len(rec_coeffs) - 1

    def munchetal_filter(im, wlevel, sigma, wname='db15'):
        # Wavelet decomposition:
        coeffs = pywt.wavedec2(im.astype(np.float32), wname, level=wlevel)
        coeffsFlt = [coeffs[0]]
        # FFT transform of horizontal frequency bands:
        for i in range(1, wlevel + 1):
            # FFT:
            fcV = np.fft.fftshift(np.fft.fft(coeffs[i][1], axis=0))
            my, mx = fcV.shape
            # Damping of vertical stripes:
            damp = 1 - np.exp(-(np.arange(-np.floor(my / 2.), -np.floor(my / 2.) + my) ** 2) / (2 * (sigma ** 2)))
            dampprime = np.kron(np.ones((1, mx)), damp.reshape((damp.shape[0], 1)))
            fcV = fcV * dampprime
            # Inverse FFT:
            fcVflt = np.real(np.fft.ifft(np.fft.ifftshift(fcV), axis=0))
            cVHDtup = (coeffs[i][0], fcVflt, coeffs[i][2])
            coeffsFlt.append(cVHDtup)

        # Get wavelet reconstruction:
        im_f = np.real(pywt.waverec2(coeffsFlt, wname))
        # Return image according to input type:
        if (im.dtype == 'uint16'):
            # Check extrema for uint16 images:
            im_f[im_f < np.iinfo(np.uint16).min] = np.iinfo(np.uint16).min
            im_f[im_f > np.iinfo(np.uint16).max] = np.iinfo(np.uint16).max
            # Return filtered image (an additional row and/or column might be present):
            return im_f[0:im.shape[0], 0:im.shape[1]].astype(np.uint16)
        else:
            return im_f[0:im.shape[0], 0:im.shape[1]]

def haar_recomp(params, rest_coeff, final_arr):
    reverse_flatten(rest_coeff, final_arr, 0)
    final_tuple = [ params]
    for i in rest_coeff:
        final_tuple.append(i)
    final_tuple = tuple(final_tuple)
    return pywt.waverec2(final_tuple, 'haar')

def w2d(img, mode='haar', level=1):
    imArray = cv2.imread(img)
    #Datatype conversions
    #convert to grayscale
    imArray = cv2.cvtColor( imArray,cv2.COLOR_RGB2GRAY )
    #convert to float
    imArray =  np.float32(imArray)
    imArray /= 255.;
    # compute coefficients
    coeffs=pywt.wavedec2(imArray, mode, level=level)

    #print len(coeffs)

    #Process Coefficients
    coeffs_H=list(coeffs[1][0])
    coeffs_H *= 0
    coeffs[1][0] = coeffs_H

    # reconstruction
    imArray_H=pywt.waverec2(coeffs, mode);
    imArray_H *= 255.;
    imArray_H =  np.uint8(imArray_H)
    #Display result
    cv2.imshow('image',imArray_H)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

    def idwt2(self):
        """
        Test pypwt for DWT reconstruction (waverec2).
        """

        W = self.W
        levels = self.levels
        # inverse DWT with pypwt
        W.forward()
        logging.info("computing Wavelets.inverse from pypwt")
        t0 = time()
        W.inverse()
        logging.info("Wavelets.inverse took %.3f ms" % elapsed_ms(t0))

        if self.do_pywt:
            # inverse DWT with pywt
            Wpy = pywt.wavedec2(self.data, self.wname, mode=per_kw, level=levels)
            logging.info("computing waverec2 from pywt")
            _ = pywt.waverec2(Wpy, self.wname, mode=per_kw)
            logging.info("pywt took %.3f ms" % elapsed_ms(t0))

        # Check reconstruction
        W_image = W.image
        maxerr = _calc_errors(self.data, W_image, "[rec]")
        self.assertTrue(maxerr < self.tol, msg="[%s] something wrong with the reconstruction (errmax = %e)" % (self.wname, maxerr))

def blend_images(base, texture, level=4, mode='sp1', base_gain=None, texture_gain=None):
    base_data = image2array(base)
    texture_data = image2array(texture)
    output_data = []

    for base_band, texture_band in zip(base_data, texture_data):
        base_band_coeffs = pywt.wavedec2(base_band, 'db2', mode, level)
        texture_band_coeffs = pywt.wavedec2(texture_band, 'db2', mode, level)

        output_band_coeffs = [base_band_coeffs[0]]
        del base_band_coeffs[0], texture_band_coeffs[0]

        for n, (base_band_details, texture_band_details) in enumerate(
            zip(base_band_coeffs, texture_band_coeffs)):
            blended_details = []
            for (base_detail, texture_detail) in zip(base_band_details, texture_band_details):
                if base_gain is not None:
                    base_detail *= base_gain
                if texture_gain is not None:
                    texture_detail *= texture_gain

                blended = numpy.where(abs(base_detail) > abs(texture_detail), base_detail, texture_detail)
                blended_details.append(blended)

            base_band_coeffs[n] = texture_band_coeffs[n] = None
            output_band_coeffs.append(blended_details)

        new_band = pywt.waverec2(output_band_coeffs, 'db2', mode)
        output_data.append(new_band)
        del new_band, base_band_coeffs, texture_band_coeffs

    del base_data, texture_data
    output_data = numpy.array(output_data)
    return array2image(output_data, base.mode)

    def _call(self, coeff):
        """Compute the discrete 1D, 2D or 3D inverse wavelet transform.

        Parameters
        ----------
        coeff : `DiscreteLpVector`

        Returns
        -------
        arr : `DiscreteLpVector`

        """
        if len(self.range.shape) == 1:
            coeff_list = array_to_pywt_coeff(coeff, self.size_list)
            x = pywt.waverec(coeff_list, self.wbasis, self.mode)
            return self.range.element(x)
        elif len(self.range.shape) == 2:
            coeff_list = array_to_pywt_coeff(coeff, self.size_list)
            x = pywt.waverec2(coeff_list, self.wbasis, self.mode)
            return self.range.element(x)
        elif len(self.range.shape) == 3:
            coeff_dict = array_to_pywt_coeff(coeff, self.size_list)
            x = wavelet_reconstruction3d(coeff_dict, self.wbasis, self.mode,
                                         self.nscales)
            return self.range.element(x)

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 waveletDenoise(u,noiseSigma):
    wavelet = pywt.Wavelet('bior6.8')
    levels  = int( np.log2(u.shape[0]) )
    waveletCoeffs = pywt.wavedec2( u, wavelet, level=levels)
    threshold=noiseSigma*np.sqrt(2*np.log2(u.size))
    NWC = [pywt.thresholding.soft(x,threshold) for x in waveletCoeffs]
    u = pywt.waverec2( NWC, wavelet)[:u.shape[0],:u.shape[1]]
    return u

def trans(imArray, mode='haar', level=1):
    coeffs = pywt.wavedec2(imArray, mode, level=level)
    coeffs_H = list(coeffs)
    coeffs_H[0] = np.zeros(coeffs_H[0].shape)
    imArray_H = pywt.waverec2(coeffs_H, mode)
    #return imArray
    #print "img1", imArray[0]
    #print "img2", imArray_H[0]
    return imArray_H

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_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 inverse(self):
     tmp = []
     tmp.append(self.a)
     for k in range(0, len(self.h)):
         detcoefs = []
         detcoefs.append(self.h[k])
         detcoefs.append(self.v[k])
         detcoefs.append(self.d[k])
         tmp.append(detcoefs)
     res = pywt.waverec2(tmp, self.wname, mode='per')
     if self.do_random_shifts is True: return _circular_shift(res, -self.shifts[0], -self.shifts[1])
     else: return res

def w2d(img, mode='haar', level=1):
    #imArray = cv2.imread(img)
    imArray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    imArray = np.float32(imArray)
    imArray /=255
    coeffs = pywt.wavedec2(imArray, mode, level=level)
    coeffs_H = list(coeffs)
    coeffs_H[0]*=0
    imArray_H = pywt.waverec2(coeffs_H, mode)
    imArray_H *= 255
    imArray_H = np.uint8(imArray_H)
    cv2.imshow("image", imArray_H)

 def rmatvec(x):
     iinf = 0
     isup = b[0].size
     yl = [x[iinf:isup].reshape(b[0].shape), ]
     for i in xrange(1, len(b)):
         tmp = list()
         for j in xrange(3):
             iinf = copy(isup)
             isup = iinf + b[i][j].size
             tmp.append(x[iinf:isup].reshape(b[i][j].shape))
         yl.append(tmp)
     return pywt.waverec2(yl, wavelet, mode=mode)[:a.shape[0], :a.shape[1]].flatten()

def prob4(image="swanlake_polluted.jpg"):
    img = imread(image, True)
    wavelet = pywt.Wavelet('haar')
    coeffs = pywt.wavedec2(img, wavelet)
    print len(coeffs)
    cleaned = pywt.waverec2(coeffs[:-1], wavelet)
    plt.subplot(121)
    plt.imshow(img,cmap='gray')
    plt.title("Original")
    plt.subplot(122)
    plt.imshow(cleaned, cmap='gray')
    plt.title("Cleaned")
    plt.show()

def test_pywt_coeff_to_array_and_array_to_pywt_coeff():
    # Verify that the helper function does indeed work as expected
    wbasis = pywt.Wavelet('db1')
    mode = 'zpd'
    nscales = 2
    n = 16
    # 1D test
    size_list = coeff_size_list((n,), nscales, wbasis, mode)
    x = np.random.rand(n)
    coeff_list = pywt.wavedec(x, wbasis, mode, nscales)
    coeff_arr = pywt_coeff_to_array(coeff_list, size_list)
    assert isinstance(coeff_arr, (np.ndarray))
    length_of_array = np.prod(size_list[0])
    length_of_array += sum(np.prod(shape) for shape in size_list[1:-1])
    assert all_equal(len(coeff_arr), length_of_array)

    coeff_list2 = array_to_pywt_coeff(coeff_arr, size_list)
    assert all_equal(coeff_list, coeff_list2)
    reconstruction = pywt.waverec(coeff_list2, wbasis, mode)
    assert all_almost_equal(reconstruction, x)

    # 2D test
    size_list = coeff_size_list((n, n), nscales, wbasis, mode)
    x = np.random.rand(n, n)
    coeff_list = pywt.wavedec2(x, wbasis, mode, nscales)
    coeff_arr = pywt_coeff_to_array(coeff_list, size_list)
    assert isinstance(coeff_arr, (np.ndarray))
    length_of_array = np.prod(size_list[0])
    length_of_array += sum(3 * np.prod(shape) for shape in size_list[1:-1])
    assert all_equal(len(coeff_arr), length_of_array)

    coeff_list2 = array_to_pywt_coeff(coeff_arr, size_list)
    assert all_equal(coeff_list, coeff_list2)
    reconstruction = pywt.waverec2(coeff_list2, wbasis, mode)
    assert all_almost_equal(reconstruction, x)

    # 3D test
    size_list = coeff_size_list((n, n, n), nscales, wbasis, mode)
    x = np.random.rand(n, n, n)
    coeff_dict = wavelet_decomposition3d(x, wbasis, mode, nscales)
    coeff_arr = pywt_coeff_to_array(coeff_dict, size_list)
    assert isinstance(coeff_arr, (np.ndarray))
    length_of_array = np.prod(size_list[0])
    length_of_array += sum(7 * np.prod(shape) for shape in size_list[1:-1])
    assert len(coeff_arr) == length_of_array

    coeff_dict2 = array_to_pywt_coeff(coeff_arr, size_list)
    reconstruction = wavelet_reconstruction3d(coeff_dict2, wbasis, mode,
                                              nscales)
    assert all_equal(coeff_dict, coeff_dict)
    assert all_almost_equal(reconstruction, x)

 def func(dframe):
     frame1, frame2 = dframe[0], dframe[1]
     frame1 = np.array(frame1)
     frame2 = np.array(frame2)
     C = pywt.wavedec2(frame1, 'db4', level=level)
     S = pywt.wavedec2(frame2, 'db4', level=level)
     tA2 = (C[0] + S[0])/2
     coeffs = fuse(tA2, C[1:], S[1:])
     fuse_img = pywt.waverec2(coeffs, 'db4')
     if frame1.dtype == np.uint16:
         fuse_img = fuse_img.clip(0,65535).astype(np.uint16)
     elif frame1.dtype == np.uint8:
         fuse_img = fuse_img.clip(0,255).astype(np.uint8)
     return np.squeeze(fuse_img)