Inverse wavelet transform: non negligible error

[megalinier]

Hello,

First of all, thank you very much for your library, it may help me a lot!

I have run some tests on a simple image (a disk of ones on a background of zero-pixels) using your library. I find a difference between the original image myImage and the result of idwt2d(dwt2d(myImage)) equal to 3e-3 while I get an error of 1e-11 using the package PyWavelets. Please find my program below. Would you have an explanation for this difference, or a suggestion about how to correct it?

(PS: I have tried using both cyclic_conv1d and cyclic_conv1d_alt, but the error remains in the two cases.)

import tfwavelets as tfw
import pywt

size_image = 256
level_decomp = 2
#%% Creation new image
a, b = 128, 128
size_image = 256
r = 64
y,x = np.ogrid[-a:size_image-a,-b:size_image-b]
mask = x*x + y*y <= r*r
new_im = np.zeros((size_image, size_image))
new_im[mask] = 255

#%% 2d dwt nd idwt using pywt
# DECOMPOSITION
im_dwt_pywt, _ = pywt.coeffs_to_array(pywt.wavedec2(new_im, 'haar',\
          mode='periodization', level=level_decomp))
plt.figure();plt.imshow(im_dwt_pywt)

# RECONSTRUCTION
_, coeff_slices = pywt.coeffs_to_array(pywt.wavedecn(np.zeros((size_image,size_image)), 'haar', \
                                                    mode='periodization', level=level_decomp))
coeffs_from_arr = pywt.array_to_coeffs(im_dwt_pywt, coeff_slices)
im_recon_pywt = pywt.waverecn(coeffs_from_arr, wavelet='haar')
plt.figure();plt.imshow(im_recon_pywt)

#%% 2d dwt nd idwt using tfwavelets 
# DECOMPOSITION
# Placeholder for input signal
tf_signal = tf.placeholder(dtype=tf.float32, shape=(size_image,size_image,1))

# Set up tf graph
output = tfw.nodes.dwt2d(tf_signal, wavelet=tfw.dwtcoeffs.haar, levels=level_decomp)

# Compute graph
with tf.Session() as sess:
    signal = sess.run(output, feed_dict={tf_signal: new_im[:,:,np.newaxis]})

plt.figure();plt.imshow(np.squeeze(signal))

# RECONSTRUCTION
# Placeholder for input signal
tf_wave = tf.placeholder(dtype=tf.float32, shape=(size_image,size_image,1))

# Set up tf graph
output = tfw.nodes.idwt2d(tf_wave, wavelet=tfw.dwtcoeffs.haar, levels=level_decomp)

# Compute graph
with tf.Session() as sess:
    recon = sess.run(output, feed_dict={tf_wave: signal})

plt.figure();plt.imshow(np.squeeze(recon))

#%% differences
plt.figure();plt.imshow(np.squeeze(recon)-new_im)
plt.figure();plt.imshow(im_recon_pywt-new_im)

#%% errors
print(np.linalg.norm(np.squeeze(recon)-new_im))
print(np.linalg.norm(im_recon_pywt-new_im))

[mathialo]

Hello! And thanks for the kind words!

The difference is 32bit vs 64bit computations. The tfwavelet library is built on 32bit floats as thats the customary datatype for commercial GPUs. When you write np.zeros(..) like you do in your code example you get the np.float64 datatype.

This allows for larger roundoff errors, meaning the comparison you make in your code example is not quite fair. If we redo your test with np.float32 as the datatype for our image, we get the same error as for tfwavelets:

import pywt
import numpy as np

size_image = 256
level_decomp = 2

# Create new image
a, b = 128, 128
size_image = 256
r = 64
y,x = np.ogrid[-a:size_image-a,-b:size_image-b]
mask = x*x + y*y <= r*r
new_im = np.zeros((size_image, size_image), dtype=np.float32)
new_im[mask] = 255

# 2d dwt and idwt using pywt
# DECOMPOSITION
im_dwt_pywt, _ = pywt.coeffs_to_array(pywt.wavedec2(new_im, 'haar',\
          mode='periodization', level=level_decomp))

# RECONSTRUCTION
_, coeff_slices = pywt.coeffs_to_array(pywt.wavedecn(np.zeros((size_image,size_image)), 'haar', \
                                                    mode='periodization', level=level_decomp))
coeffs_from_arr = pywt.array_to_coeffs(im_dwt_pywt, coeff_slices)
im_recon_pywt = pywt.waverecn(coeffs_from_arr, wavelet='haar')

# DIFF
print(np.linalg.norm(im_recon_pywt-new_im))

which when run yields:

0.003512114

[megalinier]

Thank you!