import pywt
# one‐level 2D DWT
cA, (cH, cV, cD) = pywt.dwt2(I, 'db2', mode='periodization')
# cA = LL, cH = LH, cV = HL, cD = HH

# wavelet_image_demo_levels.py

import numpy as np
import matplotlib.pyplot as plt
import pywt
from PIL import Image

# 1) Load image as grayscale float in [0,1]
img = Image.open('i1.jpg').convert('L')
I   = np.array(img, dtype=float)/255.0

# 2) Add Gaussian noise for denoising demo
sigma_noise = 0.1
I_noisy = I + sigma_noise*np.random.randn(*I.shape)
I_noisy = np.clip(I_noisy, 0, 1)

# 3) Settings
wavelet = 'db6'
levels  = 2
mode    = 'periodization'

# 4) Full 2‐level decomposition
coeffs2 = pywt.wavedec2(I_noisy, wavelet, mode=mode, level=levels)
cA2, (cH2, cV2, cD2), (cH1b, cV1b, cD1b) = coeffs2

# 5) 1‐level decomposition to get level‐1 approximation
coeffs1 = pywt.wavedec2(I_noisy, wavelet, mode=mode, level=1)
cA1, (cH1, cV1, cD1) = coeffs1

# 6) Plot everything
fig, axes = plt.subplots(levels+1, 4, figsize=(12, 3*(levels+1)))

# Row 0: Original & Noisy
axes[0,0].imshow(I,       cmap='gray'); axes[0,0].set_title('Original'); axes[0,0].axis('off')
axes[0,1].imshow(I_noisy, cmap='gray'); axes[0,1].set_title('Noisy');    axes[0,1].axis('off')
axes[0,2].axis('off'); axes[0,3].axis('off')

# Row 1: Level-2 subbands
row = axes[1]
row[0].imshow(cA2, cmap='gray'); row[0].set_title('LL₂');   row[0].axis('off')
row[1].imshow(cH2, cmap='gray'); row[1].set_title('LH₂');   row[1].axis('off')
row[2].imshow(cV2, cmap='gray'); row[2].set_title('HL₂');   row[2].axis('off')
row[3].imshow(cD2, cmap='gray'); row[3].set_title('HH₂');   row[3].axis('off')

# Row 2: Level-1 subbands
row = axes[2]
row[0].imshow(cA1, cmap='gray'); row[0].set_title('LL₁');   row[0].axis('off')
row[1].imshow(cH1, cmap='gray'); row[1].set_title('LH₁');   row[1].axis('off')
row[2].imshow(cV1, cmap='gray'); row[2].set_title('HL₁');   row[2].axis('off')
row[3].imshow(cD1, cmap='gray'); row[3].set_title('HH₁');   row[3].axis('off')

plt.tight_layout()
plt.show()

# 7) Denoising: threshold detail coefficients (universal soft)
sigma_est = np.median(np.abs(cD1)) / 0.6745
T         = sigma_est * np.sqrt(2*np.log(I.size))

# apply soft threshold to **all** detail bands at both levels
thr_coeffs2 = [cA2,
               (pywt.threshold(cH2,T,'soft'),
                pywt.threshold(cV2,T,'soft'),
                pywt.threshold(cD2,T,'soft')),
               (pywt.threshold(cH1b,T,'soft'),
                pywt.threshold(cV1b,T,'soft'),
                pywt.threshold(cD1b,T,'soft'))]

I_denoised = pywt.waverec2(thr_coeffs2, wavelet, mode=mode)
I_denoised = np.clip(I_denoised, 0, 1)

# 8) Show noisy vs. denoised
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.imshow(I_noisy, cmap='gray');    plt.title(f'Noisy (σ={sigma_noise})'); plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(I_denoised, cmap='gray'); plt.title(f'Denoised (T={T:.3f})');    plt.axis('off')
plt.tight_layout()
plt.show()

Output:

Future Work