Thresholding and Denoising

def threshold(data, value, mode: str = 'soft', substitute: float = 0):
    """
    Apply thresholding to data.
    
    Parameters:
    - data: Input array to threshold
    - value: Threshold value (scalar or array-like for per-element thresholds)
    - mode: Thresholding mode ('soft', 'hard', 'garrote', 'greater', 'less')  
    - substitute: Substitute value for thresholded elements (default: 0)
    
    Returns:
    Thresholded array of same shape as input
    """

import pywt
import numpy as np
import matplotlib.pyplot as plt

# Create test signal with noise
np.random.seed(42)
t = np.linspace(0, 1, 1000)
clean_signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise = 0.5 * np.random.randn(len(t))
noisy_signal = clean_signal + noise

print(f"Signal length: {len(noisy_signal)}")
print(f"Noise std: {np.std(noise):.3f}")

# Wavelet decomposition
coeffs = pywt.wavedec(noisy_signal, 'db8', level=6)
print(f"Decomposition levels: {len(coeffs) - 1}")

# Different thresholding modes
threshold_value = 0.1
modes = ['soft', 'hard', 'garrote', 'greater', 'less']

# Apply different thresholding modes to detail coefficients
thresholded_results = {}

for mode in modes:
    coeffs_thresh = coeffs.copy()
    
    # Apply thresholding to all detail coefficients (skip approximation)
    for i in range(1, len(coeffs_thresh)):
        coeffs_thresh[i] = pywt.threshold(coeffs_thresh[i], threshold_value, mode=mode)
    
    # Reconstruct signal
    reconstructed = pywt.waverec(coeffs_thresh, 'db8')
    thresholded_results[mode] = reconstructed

# Visualization
fig, axes = plt.subplots(3, 2, figsize=(15, 12))
axes = axes.flatten()

# Original signals
axes[0].plot(t, clean_signal, 'b-', label='Clean signal', linewidth=2)
axes[0].plot(t, noisy_signal, 'r-', alpha=0.7, label='Noisy signal')
axes[0].set_title('Original Signals')
axes[0].legend()
axes[0].grid(True)

# Thresholded results
for i, (mode, result) in enumerate(thresholded_results.items()):
    axes[i+1].plot(t, clean_signal, 'b-', alpha=0.7, label='Clean')
    axes[i+1].plot(t, result, 'g-', linewidth=2, label=f'{mode.capitalize()} threshold')
    axes[i+1].set_title(f'{mode.capitalize()} Thresholding')
    axes[i+1].legend()
    axes[i+1].grid(True)

plt.tight_layout()
plt.show()

# Compare thresholding effectiveness
def calculate_mse(signal1, signal2):
    """Calculate Mean Squared Error."""
    return np.mean((signal1 - signal2)**2)

def calculate_snr(clean, noisy):
    """Calculate Signal-to-Noise Ratio in dB."""
    signal_power = np.mean(clean**2)
    noise_power = np.mean((noisy - clean)**2)
    return 10 * np.log10(signal_power / noise_power)

print("\nThresholding Performance:")
print(f"Original SNR: {calculate_snr(clean_signal, noisy_signal):.2f} dB")
print()

for mode, result in thresholded_results.items():
    mse = calculate_mse(clean_signal, result)
    snr = calculate_snr(clean_signal, result)
    print(f"{mode.capitalize():>8}: MSE = {mse:.6f}, SNR = {snr:.2f} dB")

# Demonstrate coefficient-specific thresholding
coeffs_detailed = pywt.wavedec(noisy_signal, 'db8', level=6)

# Analyze coefficient statistics
print("\nCoefficient Statistics:")
for i, coeff in enumerate(coeffs_detailed):
    coeff_type = "Approximation" if i == 0 else f"Detail {len(coeffs_detailed)-i}"
    print(f"{coeff_type:>15}: std = {np.std(coeff):.4f}, max = {np.max(np.abs(coeff)):.4f}")

# Adaptive thresholding based on coefficient statistics
coeffs_adaptive = coeffs_detailed.copy()
for i in range(1, len(coeffs_adaptive)):  # Skip approximation
    # Use different threshold for each level
    sigma = np.std(coeffs_adaptive[i])
    adaptive_threshold = 2 * sigma  # 2-sigma rule
    coeffs_adaptive[i] = pywt.threshold(coeffs_adaptive[i], adaptive_threshold, mode='soft')

adaptive_result = pywt.waverec(coeffs_adaptive, 'db8')
adaptive_snr = calculate_snr(clean_signal, adaptive_result)
print(f"\nAdaptive thresholding SNR: {adaptive_snr:.2f} dB")

def threshold_firm(data, value_low, value_high):
    """
    Firm thresholding (intermediate between soft and hard).
    
    Parameters:
    - data: Input array to threshold
    - value_low: Lower threshold value
    - value_high: Upper threshold value (must be >= value_low)
    
    Returns:
    Firm thresholded array
    """

import pywt
import numpy as np
import matplotlib.pyplot as plt

# Create test signal with mixed noise characteristics
t = np.linspace(0, 2, 2000)
signal = (np.sin(2 * np.pi * 3 * t) +           # Low frequency component
          0.7 * np.sin(2 * np.pi * 15 * t) +    # Medium frequency component
          0.4 * np.sin(2 * np.pi * 50 * t))     # High frequency component

# Add different types of noise
gaussian_noise = 0.3 * np.random.randn(len(t))
impulse_noise = np.zeros_like(t)
impulse_indices = np.random.choice(len(t), size=50, replace=False)
impulse_noise[impulse_indices] = 2 * np.random.randn(50)

mixed_noisy = signal + gaussian_noise + impulse_noise

# Wavelet decomposition
coeffs = pywt.wavedec(mixed_noisy, 'db6', level=7)

# Compare soft, hard, and firm thresholding
threshold_low = 0.05
threshold_high = 0.15

# Soft thresholding
coeffs_soft = coeffs.copy()
for i in range(1, len(coeffs_soft)):
    coeffs_soft[i] = pywt.threshold(coeffs_soft[i], threshold_high, mode='soft')
reconstructed_soft = pywt.waverec(coeffs_soft, 'db6')

# Hard thresholding  
coeffs_hard = coeffs.copy()
for i in range(1, len(coeffs_hard)):
    coeffs_hard[i] = pywt.threshold(coeffs_hard[i], threshold_high, mode='hard')
reconstructed_hard = pywt.waverec(coeffs_hard, 'db6')

# Firm thresholding
coeffs_firm = coeffs.copy()
for i in range(1, len(coeffs_firm)):
    coeffs_firm[i] = pywt.threshold_firm(coeffs_firm[i], threshold_low, threshold_high)
reconstructed_firm = pywt.waverec(coeffs_firm, 'db6')

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

axes[0, 0].plot(t, signal, 'b-', label='Clean signal', linewidth=2)
axes[0, 0].plot(t, mixed_noisy, 'r-', alpha=0.7, label='Noisy signal')
axes[0, 0].set_title('Original vs Noisy Signal')
axes[0, 0].legend()
axes[0, 0].grid(True)

axes[0, 1].plot(t, signal, 'b-', alpha=0.7, label='Clean')
axes[0, 1].plot(t, reconstructed_soft, 'g-', linewidth=2, label='Soft threshold')
axes[0, 1].set_title('Soft Thresholding')
axes[0, 1].legend()
axes[0, 1].grid(True)

axes[1, 0].plot(t, signal, 'b-', alpha=0.7, label='Clean')
axes[1, 0].plot(t, reconstructed_hard, 'r-', linewidth=2, label='Hard threshold')
axes[1, 0].set_title('Hard Thresholding')
axes[1, 0].legend()
axes[1, 0].grid(True)

axes[1, 1].plot(t, signal, 'b-', alpha=0.7, label='Clean')
axes[1, 1].plot(t, reconstructed_firm, 'm-', linewidth=2, label='Firm threshold')
axes[1, 1].set_title('Firm Thresholding')
axes[1, 1].legend()
axes[1, 1].grid(True)

plt.tight_layout()
plt.show()

# Performance comparison
def calculate_metrics(clean, denoised):
    """Calculate denoising metrics."""
    mse = np.mean((clean - denoised)**2)
    snr = 10 * np.log10(np.mean(clean**2) / np.mean((clean - denoised)**2))
    
    # Peak Signal-to-Noise Ratio
    max_val = np.max(clean)
    psnr = 20 * np.log10(max_val / np.sqrt(mse))
    
    return mse, snr, psnr

print("Thresholding Method Comparison:")
print(f"{'Method':<15} {'MSE':<10} {'SNR (dB)':<10} {'PSNR (dB)':<10}")
print("-" * 50)

original_mse, original_snr, original_psnr = calculate_metrics(signal, mixed_noisy)
print(f"{'Original':<15} {original_mse:<10.6f} {original_snr:<10.2f} {original_psnr:<10.2f}")

soft_mse, soft_snr, soft_psnr = calculate_metrics(signal, reconstructed_soft)
print(f"{'Soft':<15} {soft_mse:<10.6f} {soft_snr:<10.2f} {soft_psnr:<10.2f}")

hard_mse, hard_snr, hard_psnr = calculate_metrics(signal, reconstructed_hard)
print(f"{'Hard':<15} {hard_mse:<10.6f} {hard_snr:<10.2f} {hard_psnr:<10.2f}")

firm_mse, firm_snr, firm_psnr = calculate_metrics(signal, reconstructed_firm)
print(f"{'Firm':<15} {firm_mse:<10.6f} {firm_snr:<10.2f} {firm_psnr:<10.2f}")

# Demonstrate threshold selection strategies
def sure_threshold_estimate(coeffs, sigma=None):
    """Estimate threshold using Stein's Unbiased Risk Estimate (SURE)."""
    if sigma is None:
        # Estimate noise standard deviation from finest detail coefficients
        sigma = np.median(np.abs(coeffs[-1])) / 0.6745
    
    n = len(coeffs[-1])
    threshold = sigma * np.sqrt(2 * np.log(n))
    return threshold

def bayes_threshold_estimate(coeffs):
    """Simple Bayesian threshold estimate."""
    # Estimate based on coefficient statistics
    all_details = np.concatenate([coeff for coeff in coeffs[1:]])
    return np.std(all_details)

# Apply different threshold selection methods
sure_thresh = sure_threshold_estimate(coeffs)
bayes_thresh = bayes_threshold_estimate(coeffs)

print(f"\nThreshold Selection:")
print(f"SURE threshold: {sure_thresh:.4f}")
print(f"Bayesian threshold: {bayes_thresh:.4f}")
print(f"Manual threshold: {threshold_high:.4f}")

# Test SURE threshold
coeffs_sure = coeffs.copy()
for i in range(1, len(coeffs_sure)):
    coeffs_sure[i] = pywt.threshold(coeffs_sure[i], sure_thresh, mode='soft')
reconstructed_sure = pywt.waverec(coeffs_sure, 'db6')

sure_mse, sure_snr, sure_psnr = calculate_metrics(signal, reconstructed_sure)
print(f"{'SURE':<15} {sure_mse:<10.6f} {sure_snr:<10.2f} {sure_psnr:<10.2f}")

import pywt
import numpy as np
import matplotlib.pyplot as plt

# Create test image
size = 256
x, y = np.mgrid[0:size, 0:size]
image = np.zeros((size, size))

# Add geometric shapes
image[64:192, 64:192] = 1.0  # Square
center_x, center_y = 128, 128
radius = 40
circle_mask = (x - center_x)**2 + (y - center_y)**2 <= radius**2
image[circle_mask] = 0.5

# Add texture and noise
texture = 0.1 * np.sin(2 * np.pi * x / 16) * np.cos(2 * np.pi * y / 16)
noise = 0.3 * np.random.randn(size, size)  
noisy_image = image + texture + noise

print(f"Image shape: {noisy_image.shape}")
print(f"Pixel value range: [{np.min(noisy_image):.3f}, {np.max(noisy_image):.3f}]")

# 2D wavelet decomposition
coeffs_2d = pywt.wavedec2(noisy_image, 'db4', level=4)
print(f"2D decomposition levels: {len(coeffs_2d) - 1}")

# Apply thresholding to 2D coefficients
def threshold_2d_coeffs(coeffs, threshold_value, mode='soft'):
    """Apply thresholding to 2D wavelet coefficients."""
    coeffs_thresh = [coeffs[0]]  # Keep approximation unchanged
    
    for level_coeffs in coeffs[1:]:
        cH, cV, cD = level_coeffs
        cH_thresh = pywt.threshold(cH, threshold_value, mode=mode)
        cV_thresh = pywt.threshold(cV, threshold_value, mode=mode)
        cD_thresh = pywt.threshold(cD, threshold_value, mode=mode)
        coeffs_thresh.append((cH_thresh, cV_thresh, cD_thresh))
    
    return coeffs_thresh

# Different threshold values
thresholds = [0.05, 0.1, 0.2, 0.3]
denoised_images = {}

for thresh in thresholds:
    coeffs_thresh = threshold_2d_coeffs(coeffs_2d, thresh, mode='soft')
    denoised = pywt.waverec2(coeffs_thresh, 'db4')
    denoised_images[thresh] = denoised

# Visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))

axes[0, 0].imshow(image, cmap='gray')
axes[0, 0].set_title('Original Clean Image')
axes[0, 0].axis('off')

axes[0, 1].imshow(noisy_image, cmap='gray')
axes[0, 1].set_title('Noisy Image')
axes[0, 1].axis('off')

# Show wavelet decomposition
cA = coeffs_2d[0]
axes[0, 2].imshow(cA, cmap='gray')
axes[0, 2].set_title('Approximation (Level 4)')
axes[0, 2].axis('off')

# Show denoised results
for i, (thresh, denoised) in enumerate(denoised_images.items()):
    if i < 3:
        axes[1, i].imshow(denoised, cmap='gray')
        axes[1, i].set_title(f'Denoised (threshold = {thresh})')
        axes[1, i].axis('off')

plt.tight_layout()
plt.show()

# Quantitative evaluation
def image_quality_metrics(original, denoised):
    """Calculate image quality metrics."""
    mse = np.mean((original - denoised)**2)
    
    # Peak Signal-to-Noise Ratio
    max_val = np.max(original)
    psnr = 20 * np.log10(max_val / np.sqrt(mse)) if mse > 0 else float('inf')
    
    # Structural Similarity Index (simplified)
    mu1, mu2 = np.mean(original), np.mean(denoised)
    sigma1, sigma2 = np.std(original), np.std(denoised)
    sigma12 = np.mean((original - mu1) * (denoised - mu2))
    
    c1, c2 = 0.01**2, 0.03**2  # Constants for stability
    ssim = ((2 * mu1 * mu2 + c1) * (2 * sigma12 + c2)) / ((mu1**2 + mu2**2 + c1) * (sigma1**2 + sigma2**2 + c2))
    
    return mse, psnr, ssim

print("\n2D Image Denoising Results:")
print(f"{'Threshold':<10} {'MSE':<10} {'PSNR (dB)':<12} {'SSIM':<10}")
print("-" * 50)

for thresh, denoised in denoised_images.items():
    mse, psnr, ssim = image_quality_metrics(image, denoised)
    print(f"{thresh:<10} {mse:<10.6f} {psnr:<12.2f} {ssim:<10.4f}")

# Adaptive 2D thresholding
def adaptive_2d_threshold(coeffs):
    """Apply adaptive thresholding based on coefficient statistics at each level."""
    coeffs_adaptive = [coeffs[0]]  # Keep approximation
    
    for level, (cH, cV, cD) in enumerate(coeffs[1:]):
        # Calculate adaptive thresholds for each orientation
        thresh_H = 3 * np.std(cH)
        thresh_V = 3 * np.std(cV) 
        thresh_D = 3 * np.std(cD)
        
        cH_thresh = pywt.threshold(cH, thresh_H, mode='soft')
        cV_thresh = pywt.threshold(cV, thresh_V, mode='soft')
        cD_thresh = pywt.threshold(cD, thresh_D, mode='soft')
        
        coeffs_adaptive.append((cH_thresh, cV_thresh, cD_thresh))
        
        print(f"Level {len(coeffs)-1-level}: thresholds H={thresh_H:.4f}, V={thresh_V:.4f}, D={thresh_D:.4f}")
    
    return coeffs_adaptive

coeffs_adaptive = adaptive_2d_threshold(coeffs_2d)
adaptive_denoised = pywt.waverec2(coeffs_adaptive, 'db4')

adaptive_mse, adaptive_psnr, adaptive_ssim = image_quality_metrics(image, adaptive_denoised)
print(f"{'Adaptive':<10} {adaptive_mse:<10.6f} {adaptive_psnr:<12.2f} {adaptive_ssim:<10.4f}")

# Show adaptive result
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.imshow(image, cmap='gray')
plt.title('Original')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.imshow(noisy_image, cmap='gray')
plt.title('Noisy')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.imshow(adaptive_denoised, cmap='gray')
plt.title('Adaptive Denoised')
plt.axis('off')

plt.tight_layout()
plt.show()

# Thresholding modes
ThresholdMode = Literal['soft', 'hard', 'garrote', 'greater', 'less']

# Threshold value types
ThresholdValue = Union[float, np.ndarray]  # Scalar or array for per-element thresholds

# Common threshold selection methods
ThresholdMethod = Literal['sure', 'bayes', 'minimax', 'manual']