Shearlet变换在图像随机噪声压制中的应用

Shearlet变换作为一种优越的多尺度几何分析工具，在图像去噪领域表现出色，特别适用于处理随机噪声。下面我将详细介绍其原理和完整的MATLAB实现方案。

Shearlet变换去噪优势

与传统方法对比

方法	优点	缺点
Shearlet变换	优越的边缘保持能力、多方向捕捉、最优稀疏表示	计算复杂度较高
小波变换	计算快速、理论基础完善	方向选择性有限、边缘保持差
傅里叶变换	全局频率分析	空间定位能力差
中值滤波	脉冲噪声抑制好	细节模糊严重

Shearlet变换去噪原理

核心思想

噪声图像 → Shearlet变换 → 系数阈值处理 → 逆Shearlet变换 → 去噪图像

噪声模型

• 加性高斯白噪声：I_noisy = I_clean + η, 其中 η ~ N(0, σ²)
• 乘性噪声：可通过对数变换转为加性噪声

MATLAB实现

1. 主函数框架

function [denoised_img, performance] = shearlet_denoise(noisy_img, varargin)
% SHEARLET_DENOISE 基于Shearlet变换的图像去噪
% 输入：
%   noisy_img - 含噪图像
%  可选参数：
%   'noise_std' - 噪声标准差估计
%   'threshold_type' - 阈值类型 ('soft', 'hard', 'bayes')
%   'scales' - 分解尺度数
%   'directional_filters' - 方向滤波器数
% 输出：
%   denoised_img - 去噪图像
%   performance - 性能指标

    % 参数解析
    p = inputParser;
    addParameter(p, 'noise_std', [], @isnumeric);
    addParameter(p, 'threshold_type', 'bayes', @ischar);
    addParameter(p, 'scales', 4, @isnumeric);
    addParameter(p, 'directional_filters', [2 3 4 4], @isnumeric);
    addParameter(p, 'display_results', true, @islogical);
    parse(p, varargin{:});
    params = p.Results;
    
    fprintf('开始Shearlet变换去噪处理...\n');
    
    % 1. 图像预处理
    preprocessed_img = preprocess_image(noisy_img);
    
    % 2. 噪声估计（如果未提供）
    if isempty(params.noise_std)
        params.noise_std = estimate_noise_std(preprocessed_img);
    end
    fprintf('估计噪声标准差: %.4f\n', params.noise_std);
    
    % 3. Shearlet变换
    shearlet_coeffs = compute_shearlet_transform(preprocessed_img, params);
    
    % 4. 系数阈值处理
    denoised_coeffs = apply_threshold(shearlet_coeffs, params);
    
    % 5. 逆Shearlet变换
    denoised_img = inverse_shearlet_transform(denoised_coeffs, params);
    
    % 6. 后处理
    denoised_img = postprocess_image(denoised_img, preprocessed_img);
    
    % 7. 性能评估
    performance = evaluate_performance(preprocessed_img, denoised_img, params);
    
    % 8. 结果可视化
    if params.display_results
        visualize_results(preprocessed_img, denoised_img, shearlet_coeffs, ...
                         denoised_coeffs, performance, params);
    end
    
    fprintf('去噪完成: PSNR=%.2fdB, SSIM=%.4f\n', ...
        performance.psnr, performance.ssim);
end

2. 图像预处理

function processed_img = preprocess_image(img)
% 图像预处理
    
    % 转换为double类型
    if ~isa(img, 'double')
        img = im2double(img);
    end
    
    % 如果是彩色图像，转换为灰度
    if ndims(img) == 3
        img = rgb2gray(img);
        fprintf('转换为灰度图像处理\n');
    end
    
    % 图像尺寸调整（优化Shearlet变换性能）
    [h, w] = size(img);
    new_size = [2^nextpow2(h), 2^nextpow2(w)];
    
    if any(new_size ~= [h, w])
        processed_img = imresize(img, new_size);
        fprintf('图像尺寸调整: %dx%d -> %dx%d\n', h, w, new_size(1), new_size(2));
    else
        processed_img = img;
    end
end

3. 噪声估计

function noise_std = estimate_noise_std(img)
% 估计图像噪声标准差
    
    % 方法1: 基于小波变换的噪声估计
    [c, s] = wavedec2(img, 1, 'db4');
    HH = detcoef2('h', c, s, 1);
    
    % 使用HH子带的稳健估计
    noise_std = median(abs(HH(:))) / 0.6745;
    
    % 方法2: 基于均匀区域的备选估计
    try
        % 寻找可能的均匀区域
        uniform_std = estimate_from_uniform_regions(img);
        if uniform_std > 0 && abs(noise_std - uniform_std) / noise_std > 0.5
            noise_std = (noise_std + uniform_std) / 2;
        end
    catch
        % 如果备选方法失败，使用小波估计
    end
    
    fprintf('噪声标准差估计: %.4f\n', noise_std);
end

function uniform_std = estimate_from_uniform_regions(img)
% 从均匀区域估计噪声
    
    % 计算局部方差
    local_var = stdfilt(img, ones(3)).^2;
    
    % 寻找低方差区域（可能是均匀区域）
    sorted_var = sort(local_var(:));
    n_lowest = round(0.1 * numel(sorted_var));
    low_var_regions = sorted_var(1:n_lowest);
    
    % 估计噪声标准差
    uniform_std = sqrt(mean(low_var_regions));
end

4. Shearlet变换实现

function shearlet_coeffs = compute_shearlet_transform(img, params)
% 计算Shearlet变换系数
    
    fprintf('计算Shearlet变换...\n');
    tic;
    
    [h, w] = size(img);
    scales = params.scales;
    directional_filters = params.directional_filters;
    
    % 初始化系数结构
    shearlet_coeffs = struct();
    shearlet_coeffs.low_pass = zeros(h, w);
    shearlet_coeffs.high_pass = cell(1, scales);
    
    % 傅里叶域处理
    img_freq = fft2(img);
    
    % 生成Shearlet滤波器
    shearlet_filters = generate_shearlet_filters(h, w, scales, directional_filters);
    
    % 应用Shearlet变换
    all_coeffs = zeros(h, w, sum(directional_filters) + 2);
    coeff_idx = 1;
    
    % 低频分量
    low_pass_coeff = ifft2(img_freq .* shearlet_filters.low_pass);
    shearlet_coeffs.low_pass = real(low_pass_coeff);
    all_coeffs(:, :, coeff_idx) = shearlet_coeffs.low_pass;
    coeff_idx = coeff_idx + 1;
    
    % 高频分量（各尺度各方向）
    for scale = 1:scales
        num_directions = directional_filters(scale);
        scale_coeffs = zeros(h, w, num_directions);
        
        for direction = 1:num_directions
            % 应用方向滤波器
            directional_filter = shearlet_filters.high_pass{scale}(:, :, direction);
            coeff = ifft2(img_freq .* directional_filter);
            scale_coeffs(:, :, direction) = real(coeff);
            
            all_coeffs(:, :, coeff_idx) = scale_coeffs(:, :, direction);
            coeff_idx = coeff_idx + 1;
        end
        
        shearlet_coeffs.high_pass{scale} = scale_coeffs;
    end
    
    shearlet_coeffs.all_coeffs = all_coeffs;
    shearlet_coeffs.filters = shearlet_filters;
    
    transform_time = toc;
    fprintf('Shearlet变换完成，耗时: %.2f秒\n', transform_time);
end

function filters = generate_shearlet_filters(h, w, scales, directional_filters)
% 生成Shearlet滤波器组
    
    % 创建频率网格
    [X, Y] = meshgrid(-w/2:w/2-1, -h/2:h/2-1);
    X = fftshift(X) / (w/2);
    Y = fftshift(Y) / (h/2);
    
    R = sqrt(X.^2 + Y.^2);  % 径向频率
    Theta = atan2(Y, X);    % 角度
    
    filters = struct();
    filters.low_pass = zeros(h, w);
    filters.high_pass = cell(1, scales);
    
    % 生成Meyer小波尺度函数（低频）
    filters.low_pass = meyer_scale_function(R);
    
    % 生成各尺度方向滤波器
    for scale = 1:scales
        num_directions = directional_filters(scale);
        scale_filters = zeros(h, w, num_directions);
        
        % 尺度相关的径向窗函数
        radial_window = meyer_wavelet_function(2^(scale-1) * R);
        
        for direction = 1:num_directions
            % 方向角度
            theta_k = (direction - 1) * 2 * pi / num_directions;
            
            % 方向窗函数
            angular_window = directional_meyer_function(...
                num_directions * (Theta - theta_k));
            
            % 组合得到Shearlet滤波器
            shearlet_filter = radial_window .* angular_window;
            
            % 归一化
            scale_filters(:, :, direction) = shearlet_filter / norm(shearlet_filter(:));
        end
        
        filters.high_pass{scale} = scale_filters;
    end
end

function phi = meyer_scale_function(r)
% Meyer尺度函数
    phi = zeros(size(r));
    
    % 支撑在 [0, 2/3]
    idx1 = (r >= 0) & (r <= 2/3);
    phi(idx1) = 1;
    
    % 过渡区域 [2/3, 1]
    idx2 = (r > 2/3) & (r <= 1);
    v = r(idx2) - 2/3;
    phi(idx2) = cos(pi/2 * meyer_aux_function(3*v - 1));
end

function psi = meyer_wavelet_function(r)
% Meyer小波函数
    psi = zeros(size(r));
    
    % 支撑在 [1/3, 4/3]
    idx1 = (r >= 1/3) & (r <= 2/3);
    v1 = r(idx1) - 1/3;
    psi(idx1) = sin(pi/2 * meyer_aux_function(3*v1 - 1));
    
    idx2 = (r > 2/3) & (r <= 4/3);
    v2 = r(idx2) - 2/3;
    psi(idx2) = cos(pi/2 * meyer_aux_function(3*v2/2 - 1));
end

function w = directional_meyer_function(theta)
% 方向Meyer函数
    w = zeros(size(theta));
    
    % 支撑在 [-1, 1]
    idx = (abs(theta) <= 1);
    w(idx) = cos(pi/2 * meyer_aux_function(abs(theta(idx)) - 1/2));
end

function v = meyer_aux_function(t)
% Meyer辅助函数
    v = zeros(size(t));
    
    idx = (t >= 0) & (t <= 1);
    v(idx) = t(idx).^4 .* (35 - 84*t(idx) + 70*t(idx).^2 - 20*t(idx).^3);
end

5. 阈值处理核心算法

function denoised_coeffs = apply_threshold(shearlet_coeffs, params)
% 对Shearlet系数应用阈值处理
    
    fprintf('应用%s阈值处理...\n', params.threshold_type);
    
    denoised_coeffs = shearlet_coeffs;
    noise_std = params.noise_std;
    
    % 保留低频系数（通常包含主要图像信息）
    % 主要处理高频系数
    
    total_coeffs = 0;
    thresholded_coeffs = 0;
    
    for scale = 1:params.scales
        scale_coeffs = shearlet_coeffs.high_pass{scale};
        [h, w, num_directions] = size(scale_coeffs);
        
        for direction = 1:num_directions
            coeffs = scale_coeffs(:, :, direction);
            
            % 计算尺度相关的阈值
            threshold = calculate_threshold(coeffs, noise_std, scale, ...
                                          params.threshold_type);
            
            % 应用阈值
            thresholded_coeffs = apply_single_threshold(coeffs, threshold, ...
                                                      params.threshold_type);
            
            denoised_coeffs.high_pass{scale}(:, :, direction) = thresholded_coeffs;
            
            total_coeffs = total_coeffs + numel(coeffs);
        end
    end
    
    % 更新所有系数
    denoised_coeffs.all_coeffs(:, :, 2:end) = ...
        reconstruct_all_highpass_coeffs(denoised_coeffs);
    
    fprintf('阈值处理完成\n');
end

function threshold = calculate_threshold(coeffs, noise_std, scale, threshold_type)
% 计算阈值
    
    switch threshold_type
        case 'universal'
            % 通用阈值
            n = numel(coeffs);
            threshold = noise_std * sqrt(2 * log(n));
            
        case 'sure'
            % SURE阈值（Stein无偏风险估计）
            threshold = sure_threshold(coeffs, noise_std);
            
        case 'bayes'
            % BayesShrink阈值
            sigma_x = max(std(coeffs(:))^2 - noise_std^2, 0);
            if sigma_x > 0
                threshold = noise_std^2 / sqrt(sigma_x);
            else
                threshold = max(abs(coeffs(:)));
            end
            
        case 'scale_adaptive'
            % 尺度自适应阈值
            base_threshold = noise_std * sqrt(2 * log(numel(coeffs)));
            threshold = base_threshold / (1.5^(scale-1));
            
        otherwise
            error('未知的阈值类型: %s', threshold_type);
    end
end

function thresholded_coeffs = apply_single_threshold(coeffs, threshold, threshold_type)
% 应用单一阈值
    
    switch threshold_type
        case 'hard'
            % 硬阈值
            thresholded_coeffs = coeffs .* (abs(coeffs) > threshold);
            
        case 'soft'
            % 软阈值
            signs = sign(coeffs);
            thresholded_coeffs = signs .* max(abs(coeffs) - threshold, 0);
            
        case 'bayes'
            % Bayes软阈值
            signs = sign(coeffs);
            thresholded_coeffs = signs .* max(abs(coeffs) - threshold, 0);
            
        otherwise
            % 默认软阈值
            signs = sign(coeffs);
            thresholded_coeffs = signs .* max(abs(coeffs) - threshold, 0);
    end
end

function sure_thresh = sure_threshold(coeffs, noise_std)
% SURE阈值计算
    n = numel(coeffs);
    coeffs_sorted = sort(abs(coeffs(:)));
    
    risks = zeros(n, 1);
    for i = 1:n
        threshold = coeffs_sorted(i);
        % SURE风险计算
        num_above = sum(coeffs_sorted >= threshold);
        sum_squares = sum(min(coeffs_sorted, threshold).^2);
        risks(i) = (n - 2*num_above + sum_squares) / n;
    end
    
    [~, min_idx] = min(risks);
    sure_thresh = coeffs_sorted(min_idx);
end

6. 逆变换与后处理

function denoised_img = inverse_shearlet_transform(denoised_coeffs, params)
% Shearlet逆变换
    
    fprintf('执行Shearlet逆变换...\n');
    tic;
    
    [h, w] = size(denoised_coeffs.low_pass);
    reconstructed_freq = zeros(h, w);
    
    % 添加低频分量
    low_pass_freq = fft2(denoised_coeffs.low_pass);
    reconstructed_freq = reconstructed_freq + low_pass_freq .* ...
                        denoised_coeffs.filters.low_pass;
    
    % 添加各尺度方向的高频分量
    for scale = 1:params.scales
        scale_coeffs = denoised_coeffs.high_pass{scale};
        num_directions = size(scale_coeffs, 3);
        
        for direction = 1:num_directions
            coeff_freq = fft2(scale_coeffs(:, :, direction));
            directional_filter = denoised_coeffs.filters.high_pass{scale}(:, :, direction);
            
            reconstructed_freq = reconstructed_freq + coeff_freq .* directional_filter;
        end
    end
    
    % 逆傅里叶变换
    denoised_img = real(ifft2(reconstructed_freq));
    
    inverse_time = toc;
    fprintf('逆变换完成，耗时: %.2f秒\n', inverse_time);
end

function final_img = postprocess_image(denoised_img, original_img)
% 图像后处理
    
    % 裁剪到原始尺寸（如果之前调整过）
    [h_orig, w_orig] = size(original_img);
    [h_denoised, w_denoised] = size(denoised_img);
    
    if h_denoised ~= h_orig || w_denoised ~= w_orig
        denoised_img = imresize(denoised_img, [h_orig, w_orig]);
    end
    
    % 幅度裁剪到[0,1]
    final_img = max(0, min(1, denoised_img));
    
    % 保持与原始图像相同的动态范围
    final_img = (final_img - min(final_img(:))) / ...
                (max(final_img(:)) - min(final_img(:)));
end

7. 性能评估

function performance = evaluate_performance(original, denoised, params)
% 评估去噪性能
    
    performance = struct();
    
    % PSNR
    mse = mean((original(:) - denoised(:)).^2);
    performance.psnr = 10 * log10(1 / mse);
    
    % SSIM
    performance.ssim = ssim(denoised, original);
    
    % 边缘保持指数
    performance.epi = edge_preservation_index(original, denoised);
    
    % 特征相似性
    performance.fsim = feature_similarity_index(original, denoised);
    
    % 运行时间（如果需要）
    performance.noise_std_estimated = params.noise_std;
    
    fprintf('性能指标:\n');
    fprintf('  PSNR: %.2f dB\n', performance.psnr);
    fprintf('  SSIM: %.4f\n', performance.ssim);
    fprintf('  边缘保持指数: %.4f\n', performance.epi);
    fprintf('  特征相似性: %.4f\n', performance.fsim);
end

function epi = edge_preservation_index(original, denoised)
% 边缘保持指数
    [gx_orig, gy_orig] = gradient(original);
    [gx_den, gy_den] = gradient(denoised);
    
    grad_orig = sqrt(gx_orig.^2 + gy_orig.^2);
    grad_den = sqrt(gx_den.^2 + gy_den.^2);
    
    epi = sum(grad_orig(:) .* grad_den(:)) / ...
          (sqrt(sum(grad_orig(:).^2)) * sqrt(sum(grad_den(:).^2)));
end

function fsim = feature_similarity_index(img1, img2)
% 特征相似性指数
    % 使用相位一致性作为特征
    pc1 = phase_congruency(img1);
    pc2 = phase_congruency(img2);
    
    C = 0.001; % 常数避免除零
    fsim = (2 * pc1 .* pc2 + C) ./ (pc1.^2 + pc2.^2 + C);
    fsim = mean(fsim(:));
end

8. 结果可视化

function visualize_results(original, denoised, orig_coeffs, denoised_coeffs, ...
                          performance, params)
% 可视化去噪结果和系数
    
    figure('Position', [100, 100, 1400, 900]);
    
    % 1. 原始图像和去噪图像对比
    subplot(3, 4, 1);
    imshow(original, []);
    title('原始图像');
    
    subplot(3, 4, 2);
    imshow(denoised, []);
    title(sprintf('去噪图像\nPSNR: %.2fdB, SSIM: %.4f', ...
        performance.psnr, performance.ssim));
    
    % 2. 噪声图像（模拟）
    subplot(3, 4, 3);
    noisy_img = original + params.noise_std * randn(size(original));
    imshow(noisy_img, []);
    title(sprintf('含噪图像\n噪声标准差: %.3f', params.noise_std));
    
    % 3. 残差图像
    subplot(3, 4, 4);
    residual = original - denoised;
    imshow(residual, []);
    title('残差图像');
    colorbar;
    
    % 4. 系数分布对比（第一尺度）
    subplot(3, 4, 5);
    if params.scales >= 1
        orig_scale_coeffs = orig_coeffs.high_pass{1}(:, :, 1);
        denoised_scale_coeffs = denoised_coeffs.high_pass{1}(:, :, 1);
        
        histogram(orig_scale_coeffs(:), 50, 'Normalization', 'pdf');
        hold on;
        histogram(denoised_scale_coeffs(:), 50, 'Normalization', 'pdf');
        legend('原始系数', '去噪系数');
        title('第一尺度系数分布');
        grid on;
    end
    
    % 5. 系数能量分布
    subplot(3, 4, 6);
    scale_energy = zeros(params.scales, 1);
    for scale = 1:params.scales
        scale_coeffs = denoised_coeffs.high_pass{scale};
        scale_energy(scale) = sum(scale_coeffs(:).^2);
    end
    bar(scale_energy);
    title('各尺度系数能量');
    xlabel('尺度'); ylabel('能量');
    grid on;
    
    % 6. 边缘保持可视化
    subplot(3, 4, 7);
    orig_edges = edge(original, 'canny');
    denoised_edges = edge(denoised, 'canny');
    edge_overlap = imfuse(orig_edges, denoised_edges);
    imshow(edge_overlap);
    title(sprintf('边缘保持\nEPI: %.4f', performance.epi));
    
    % 7. 局部细节对比
    subplot(3, 4, 8);
    % 选择图像中心区域
    [h, w] = size(original);
    roi = original(round(h/4):round(3*h/4), round(w/4):round(3*w/4));
    roi_denoised = denoised(round(h/4):round(3*h/4), round(w/4):round(3*w/4));
    imshowpair(roi, roi_denoised, 'montage');
    title('局部细节对比');
    
    % 8. 频率谱分析
    subplot(3, 4, 9);
    orig_spectrum = abs(fftshift(fft2(original)));
    denoised_spectrum = abs(fftshift(fft2(denoised)));
    
    imagesc(log1p(orig_spectrum));
    title('原始图像频谱');
    colorbar;
    
    subplot(3, 4, 10);
    imagesc(log1p(denoised_spectrum));
    title('去噪图像频谱');
    colorbar;
    
    % 9. 性能指标汇总
    subplot(3, 4, 11);
    metrics = [performance.psnr, performance.ssim * 100, ...
               performance.epi * 100, performance.fsim * 100];
    metric_names = {'PSNR', 'SSIM', 'EPI', 'FSIM'};
    bar(metrics);
    set(gca, 'XTickLabel', metric_names);
    title('性能指标汇总');
    ylabel('得分');
    grid on;
    
    % 10. 参数信息
    subplot(3, 4, 12);
    text(0.1, 0.7, sprintf('参数设置:\n尺度数: %d\n方向数: %s\n阈值类型: %s\n噪声估计: %.4f', ...
        params.scales, mat2str(params.directional_filters), ...
        params.threshold_type, params.noise_std), 'FontSize', 10);
    text(0.1, 0.3, sprintf('处理信息:\n图像尺寸: %dx%d\n方法: Shearlet变换', ...
        size(original,1), size(original,2)), 'FontSize', 10);
    axis off;
end

📈 应用示例和测试

示例1：标准测试图像去噪

function demo_shearlet_denoise()
% Shearlet去噪演示
    
    % 读取测试图像
    img = im2double(imread('cameraman.tif'));
    
    % 添加高斯噪声
    noise_std = 0.1;
    noisy_img = img + noise_std * randn(size(img));
    
    % 不同阈值方法比较
    methods = {'soft', 'hard', 'bayes', 'scale_adaptive'};
    results = cell(length(methods), 1);
    
    figure('Position', [50, 50, 1200, 800]);
    
    for i = 1:length(methods)
        % Shearlet去噪
        [denoised, perf] = shearlet_denoise(noisy_img, ...
            'noise_std', noise_std, ...
            'threshold_type', methods{i}, ...
            'display_results', false);
        
        results{i} = struct('image', denoised, 'performance', perf);
        
        % 显示结果
        subplot(2, 3, i);
        imshow(denoised, []);
        title(sprintf('%s阈值\nPSNR: %.2fdB', methods{i}, perf.psnr));
    end
    
    % 显示原始和噪声图像
    subplot(2, 3, 5);
    imshow(img, []);
    title('原始图像');
    
    subplot(2, 3, 6);
    imshow(noisy_img, []);
    title(sprintf('噪声图像\n噪声标准差: %.3f', noise_std));
    
    % 性能比较表格
    fprintf('\n方法比较:\n');
    fprintf('%-15s %-8s %-8s %-8s\n', '方法', 'PSNR', 'SSIM', 'EPI');
    fprintf('%-15s %-8s %-8s %-8s\n', '------', '----', '----', '---');
    for i = 1:length(methods)
        p = results{i}.performance;
        fprintf('%-15s %-8.2f %-8.4f %-8.4f\n', ...
            methods{i}, p.psnr, p.ssim, p.epi);
    end
end

示例2：不同噪声水平测试

function noise_level_analysis()
% 不同噪声水平下的性能分析
    
    img = im2double(imread('lena.png'));
    if ndims(img) == 3
        img = rgb2gray(img);
    end
    
    noise_levels = [0.05, 0.1, 0.15, 0.2, 0.25];
    psnr_results = zeros(length(noise_levels), 3);
    ssim_results = zeros(length(noise_levels), 3);
    
    methods = {'shearlet', 'wavelet', 'bm3d'};
    
    for i = 1:length(noise_levels)
        noise_std = noise_levels(i);
        noisy_img = img + noise_std * randn(size(img));
        
        % Shearlet方法
        [denoised_shearlet, perf_shearlet] = shearlet_denoise(noisy_img, ...
            'noise_std', noise_std, 'display_results', false);
        
        % 小波方法比较
        denoised_wavelet = wavelet_denoise(noisy_img, noise_std);
        
        % BM3D方法（如果有工具箱）
        try
            denoised_bm3d = BM3D(noisy_img, noise_std);
        catch
            denoised_bm3d = denoised_wavelet; % 备选
        end
        
        % 计算性能指标
        psnr_results(i, 1) = perf_shearlet.psnr;
        psnr_results(i, 2) = 10 * log10(1 / mean((img(:) - denoised_wavelet(:)).^2));
        psnr_results(i, 3) = 10 * log10(1 / mean((img(:) - denoised_bm3d(:)).^2));
        
        ssim_results(i, 1) = perf_shearlet.ssim;
        ssim_results(i, 2) = ssim(denoised_wavelet, img);
        ssim_results(i, 3) = ssim(denoised_bm3d, img);
    end
    
    % 绘制性能曲线
    figure('Position', [100, 100, 1000, 400]);
    
    subplot(1, 2, 1);
    plot(noise_levels, psnr_results, 'o-', 'LineWidth', 2);
    legend('Shearlet', '小波', 'BM3D', 'Location', 'best');
    xlabel('噪声标准差');
    ylabel('PSNR (dB)');
    title('不同噪声水平的PSNR比较');
    grid on;
    
    subplot(1, 2, 2);
    plot(noise_levels, ssim_results, 's-', 'LineWidth', 2);
    legend('Shearlet', '小波', 'BM3D', 'Location', 'best');
    xlabel('噪声标准差');
    ylabel('SSIM');
    title('不同噪声水平的SSIM比较');
    grid on;
end

参考代码 shearlet变换的应用 www.youwenfan.com/contentcsk/79586.html

实用技巧和注意事项

1. 参数选择指南

function optimal_params = optimize_parameters(img_type, noise_level)
% 根据图像类型和噪声水平推荐参数
    
    optimal_params = struct();
    
    switch img_type
        case 'natural'
            optimal_params.scales = 4;
            optimal_params.directional_filters = [2, 3, 4, 4];
            optimal_params.threshold_type = 'bayes';
            
        case 'texture_rich'
            optimal_params.scales = 5;
            optimal_params.directional_filters = [3, 4, 5, 5, 5];
            optimal_params.threshold_type = 'scale_adaptive';
            
        case 'medical'
            optimal_params.scales = 3;
            optimal_params.directional_filters = [2, 3, 3];
            optimal_params.threshold_type = 'soft';
            
        otherwise
            optimal_params.scales = 4;
            optimal_params.directional_filters = [2, 3, 4, 4];
            optimal_params.threshold_type = 'bayes';
    end
    
    % 根据噪声水平调整
    if noise_level > 0.2
        optimal_params.threshold_type = 'hard'; % 强噪声用硬阈值
    end
end

2. 计算优化建议

function fast_shearlet_denoise(noisy_img, noise_std)
% 快速Shearlet去噪（牺牲一些精度换取速度）
    
    % 使用较小的尺度数和方向数
    params.scales = 3;
    params.directional_filters = [2, 2, 3];
    params.noise_std = noise_std;
    params.threshold_type = 'universal'; % 计算更简单
    params.display_results = true;
    
    % 下采样处理（对大图像）
    if min(size(noisy_img)) > 512
        scale_factor = 512 / min(size(noisy_img));
        small_img = imresize(noisy_img, scale_factor);
        
        [denoised_small, perf] = shearlet_denoise(small_img, params);
        
        % 上采样回原尺寸
        denoised_img = imresize(denoised_small, size(noisy_img));
    else
        [denoised_img, perf] = shearlet_denoise(noisy_img, params);
    end
    
    return denoised_img, perf;
end

扩展应用

1. 彩色图像去噪

function color_denoise(rgb_img, noise_std)
% 彩色图像Shearlet去噪
    
    % 转换到YUV色彩空间
    yuv_img = rgb2ycbcr(rgb_img);
    
    % 仅对亮度分量去噪（保持色彩信息）
    Y = yuv_img(:, :, 1);
    U = yuv_img(:, :, 2);
    V = yuv_img(:, :, 3);
    
    % 对Y分量进行Shearlet去噪
    Y_denoised = shearlet_denoise(Y, 'noise_std', noise_std, ...
                                 'display_results', false);
    
    % 重组图像
    yuv_denoised = cat(3, Y_denoised, U, V);
    rgb_denoised = ycbcr2rgb(yuv_denoised);
    
    return rgb_denoised;
end

2. 视频去噪框架

function video_denoise(video_path, output_path)
% 视频序列Shearlet去噪
    
    video_reader = VideoReader(video_path);
    video_writer = VideoWriter(output_path, 'MPEG-4');
    open(video_writer);
    
    frame_count = 0;
    while hasFrame(video_reader)
        frame = readFrame(video_reader);
        frame_double = im2double(frame);
        
        % 估计噪声（使用第一帧或连续估计）
        if frame_count == 0
            noise_std = estimate_noise_std(rgb2gray(frame_double));
        end
        
        % 彩色去噪
        denoised_frame = color_denoise(frame_double, noise_std);
        
        writeVideo(video_writer, denoised_frame);
        frame_count = frame_count + 1;
        
        if mod(frame_count, 30) == 0
            fprintf('已处理 %d 帧\n', frame_count);
        end
    end
    
    close(video_writer);
    fprintf('视频去噪完成: %d 帧\n', frame_count);
end

这个完整的Shearlet变换去噪实现提供了从理论到实践的全面解决方案。你可以根据具体图像特性和噪声水平调整参数，获得最佳的去噪效果。