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?
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?
1 Answer
Reset to default 0There is a function skimage.restoration.calibrate_denoiser that you can use to calibrate any denoising function. There is a short tutorial on how to use it here and a more in-depth one here.
noise_std*3. Please correct me if I'm wrong ;) Actually, I thought the optimal weight should be similar if the optimal value is within the weight range. – zxdawn Commented Jan 18 at 20:46