### 使用Python实现小波阈值去噪 #### 导入必要的库 为了执行小波去噪操作，需要导入一些特定的库来支持图像读取、显示以及小波变换等功能。 ```python import pywt import numpy as np from PIL import Image import matplotlib.pyplot as plt ``` #### 加载并预处理图像数据 加载待处理的图片文件，并将其转换成灰度模式以便后续处理。这里假设有一个名为`noisy_image.png`的含噪声输入图像。 ```python def load_and_preprocess(image_path): img = Image.open(image_path).convert('L') # 转换成灰度图 img_array = np.array(img, dtype=np.float64) / 255.0 # 归一化到[0,1]区间 return img_array ``` #### 执行离散小波变换(DWT) 应用多级DWT分解原图像，在此过程中会产生一系列低频子带(LL)和高频子带(LH, HL, HH)，其中LL代表近似系数而其他三个则对应细节部分。 ```python def apply_dwt(image, wavelet='db4', level=3): coeffs = pywt.wavedec2(data=image, wavelet=wavelet, level=level) cA, (cH, cV, cD) = coeffs[0], tuple(coeffs[-level:]) return cA, (cH, cV, cD), coeffs ``` #### 阈值量化(Thresholding) 对于得到的小波域表示形式，可以通过设定合适的阈值λ对各层细节成分实施硬/软阈值处理以达到降噪目的；本例采用简单的全局统一阈值策略。 ```python def threshold_coeffs(detail_coeffs, threshold_type='soft'): thresh = np.std(detail_coeffs) * np.sqrt(2*np.log(len(detail_coeffs))) if threshold_type == 'hard': detail_coeffs[np.abs(detail_coeffs) < thresh] = 0 elif threshold_type == 'soft': shrinked = np.sign(detail_coeffs) * np.maximum(np.abs(detail_coeffs)-thresh, 0) detail_coeffs[:] = shrinked return detail_coeffs ``` #### 进行逆向离散小波变换(IDWT) 完成上述步骤之后再利用IDWT重构最终净化后的版本。 ```python def reconstruct_from_dwt(cA, details, original_shape, wavelet='db4'): recon_coeffs = [cA] + list(details)[::-1] reconstructed_img = pywt.waverec2(recon_coeffs, wavelet)[:original_shape[0], :original_shape[1]] normalized_reconstructed_img = ((reconstructed_img - reconstructed_img.min()) / (reconstructed_img.max() - reconstructed_img.min())) return normalized_reconstructed_img ``` #### 完整过程封装函数 最后可以将整个流程打包成一个易于调用的功能接口供外部程序直接使用。 ```python def denoise_wavelets(input_image_path, output_image_path=None, wavelet='db4', levels=3, threshold_type='soft'): image_data = load_and_preprocess(input_image_path) _, details, all_coeffs = apply_dwt(image_data, wavelet, levels) for i in range(levels): current_level_details = all_coeffs[-i-1][:] updated_current_level_details = [] for d in current_level_details: updated_detail_coefficients = threshold_coeffs(d.flatten(), threshold_type) updated_current_level_details.append(updated_detail_coefficients.reshape(d.shape)) all_coeffs[-i-1] = tuple(updated_current_level_details) final_cA = all_coeffs[0] cleaned_image = reconstruct_from_dwt(final_cA, all_coeffs[1:], image_data.shape) fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 4)) ax[0].imshow(image_data, cmap="gray") ax[0].set_title("Original Noisy Image") ax[1].imshow(cleaned_image, cmap="gray") ax[1].set_title(f"Denoised ({threshold_type} Threshold)") for a in ax.flat: a.axis('off') plt.tight_layout() if output_image_path is not None: plt.savefig(output_image_path) else: plt.show() denoise_wavelets('path_to_noisy_image.png') ```