I have 2d data which has background noise and assembled high values. I'm trying to apply the TV filter to denoise the data. Is there a suitable method to avoid over-denoising the data? I have tried to check the MSE value like this: import numpy as np from skimage.restoration import denoise_tv_chambolle import matplotlib.pyplot as plt from sklearn.metrics import mean_squared_error from scipy import stats def create_mask_from_quantiles(image, lower_quantile=0.1, upper_quantile=0.9): """ Create a mask based on quantile values to exclude extreme values. Parameters: ----------- image : 2D numpy array Input image lower_quantile : float Lower quantile threshold (0-1) upper_quantile : float Upper quantile threshold (0-1) Returns: -------- mask : 2D boolean array True for values within quantile range """ lower_thresh = np.quantile(image, lower_quantile) upper_thresh = np.quantile(image, upper_quantile) return (image >= lower_thresh) & (image <= upper_thresh) def tune_tv_chambolle(noisy_image, weight_range=None, n_weights=10, lower_quantile=0.1, upper_quantile=0.9, multiplier=2, plot=True): """ Find suitable weight parameter for TV Chambolle denoising with quantile masking. Parameters: ----------- noisy_image : 2D numpy array Input image with noise weight_range : tuple, optional Range of weights to test (min, max) n_weights : int Number of weights to test lower_quantile : float Lower quantile for masking upper_quantile : float Upper quantile for masking multiplier: float STD*multiplier for weight_max plot : bool Whether to plot the results """ # Create mask for background noise evaluation mask = create_mask_from_quantiles(noisy_image, lower_quantile, upper_quantile) if weight_range is None: # Estimate initial range based on masked image statistics masked_image = noisy_image[mask] noise_std = np.std(masked_image) weight_range = (noise_std/10, noise_std*multiplier) print('weight_range: ', weight_range) weights = np.linspace(weight_range[0], weight_range[1], n_weights) denoised_images = [] metrics = { 'total_variation': [], 'mse_with_original': [], 'entropy': [], 'masked_mse': [] } # Calculate baseline metrics baseline_tv = np.sum(np.abs(np.diff(noisy_image, axis=0))) + \ np.sum(np.abs(np.diff(noisy_image, axis=1))) # Test different weights for weight in weights: denoised = denoise_tv_chambolle(noisy_image, weight=weight, channel_axis=None) denoised_images.append(denoised) # Calculate metrics for masked region tv = np.sum(np.abs(np.diff(denoised, axis=0))) + \ np.sum(np.abs(np.diff(denoised, axis=1))) metrics['total_variation'].append(tv) metrics['mse_with_original'].append( mean_squared_error(noisy_image, denoised)) metrics['masked_mse'].append( mean_squared_error(noisy_image[mask], denoised[mask])) metrics['entropy'].append(stats.entropy(denoised[mask].flatten())) # Find optimal weight using the elbow method on masked MSE mse_differences = np.diff(metrics['masked_mse']) elbow_idx = np.argmax(np.abs(np.diff(mse_differences))) + 1 optimal_weight = weights[elbow_idx] if plot: fig, axes = plt.subplots(2, 3, figsize=(18, 12)) # Original image im0 = axes[0,0].imshow(noisy_image, cmap='viridis') axes[0,0].set_title('Original Noisy Image') plt.colorbar(im0, ax=axes[0,0]) # Mask visualization im1 = axes[0,1].imshow(mask, cmap='gray_r') axes[0,1].set_title('Quantile Mask') plt.colorbar(im1, ax=axes[0,1]) # Denoised image optimal_denoised = denoised_images[elbow_idx] im2 = axes[0,2].imshow(optimal_denoised, cmap='viridis') axes[0,2].set_title(f'Denoised (weight={optimal_weight:.4f})') plt.colorbar(im2, ax=axes[0,2]) # Metrics plots axes[1,0].plot(weights, metrics['total_variation'], 'b-', label='Total Variation') axes[1,0].axvline(x=optimal_weight, color='r', linestyle='--', label='Optimal Weight') axes[1,0].set_xlabel('Weight') axes[1,0].set_ylabel('Total Variation') axes[1,0].legend() axes[1,1].plot(weights, metrics['masked_mse'], 'g-', label='Masked MSE') axes[1,1].set_xlabel('Weight') axes[1,1].set_ylabel('MSE (masked region)') axes[1,1].legend() axes[1,2].plot(weights, metrics['entropy'], 'm-', label='Entropy (masked region)') axes[1,2].set_xlabel('Weight') axes[1,2].set_ylabel('Entropy') axes[1,2].legend() plt.suptitle(f'multiplier={multiplier}') plt.tight_layout() plt.show() return optimal_weight, { 'weights': weights, 'metrics': metrics, 'denoised_images': denoised_images, 'mask': mask } def generate_test_data(size=100, noise_level=0.1, wind_direction=45, wind_speed=1.0): """ Generate synthetic test data with a realistic gas plume Parameters: ----------- size : int Size of the output image (size x size) noise_level : float Standard deviation of the background noise wind_direction : float Wind direction in degrees (0 is East, 90 is North) wind_speed : float Relative wind speed affecting plume spread Returns: -------- test_image : 2D numpy array Normalized image containing background noise and plume """ # Create background noise background = np.random.normal(0, noise_level, (size, size)) # Create coordinate grid x, y = np.meshgrid(np.linspace(-2, 2, size), np.linspace(-2, 2, size)) # Convert wind direction to radians and calculate rotated coordinates theta = np.radians(wind_direction) x_rot = x * np.cos(theta) + y * np.sin(theta) y_rot = -x * np.sin(theta) + y * np.cos(theta) # Create elongated plume shape # Higher wind_speed creates more elongated plume sigma_x = 0.8 * wind_speed # Spread in wind direction sigma_y = 0.2 # Spread perpendicular to wind direction # Add turbulent diffusion effect turbulence = np.random.normal(0, 0.1, (size, size)) # Calculate plume concentration with gaussian dispersion model plume = np.exp(-(x_rot**2/(2*sigma_x**2) + y_rot**2/(2*sigma_y**2))) # Add some random variations to make it more realistic plume = plume * (1 + 0.2 * turbulence) # Add source point with higher concentration source_x = size // 4 source_y = size // 2 plume[source_y-2:source_y+2, source_x-2:source_x+2] = 1.0 # Combine background and plume test_image = background + plume # Add some patchiness to the plume patchiness = np.random.normal(0, 0.05, (size, size)) test_image = test_image * (1 + patchiness * (plume > 0.1)) # Normalize to 0-1 range test_image = (test_image - test_image.min()) / (test_image.max() - test_image.min()) return test_image # Generate test data with different wind conditions test_image = generate_test_data( size=100, noise_level=0.1, wind_direction=45, # 45 degree wind direction wind_speed=1.5 # Moderate wind speed ) multiplier = 2 # Find optimal weight and denoise optimal_weight, results = tune_tv_chambolle( test_image, lower_quantile=0.2, upper_quantile=0.8, n_weights=15, multiplier=multiplier, ) # Apply final denoising with optimal weight final_denoised = denoise_tv_chambolle(test_image, weight=optimal_weight, channel_axis=None, ) # Compare original vs denoised for masked region mask = results['mask'] original_std = np.std(test_image[mask]) denoised_std = np.std(final_denoised[mask]) noise_reduction = (1 - denoised_std/original_std) * 100 print(f"Optimal weight: {optimal_weight:.4f}") print(f"Noise reduction in masked region: {noise_reduction:.1f}%") Here're some results using different weight_max of TV filter weight_max = noise_std*2 Optimal weight: 0.0572 Noise reduction in masked region: 20.8% weight_max = noise_std*3 Optimal weight: 0.0259 Noise reduction in masked region: 24.9% weight_max = noise_std*4 Optimal weight: 0.0436 Noise reduction in masked region: 25.2% Question It seems noise_std*3 gets the wrong optimal weight according to the image. How to set a correct weight range to get the optimalized TV filter? Continue reading...
