如何编码实现不同小波分解？图像去噪整合分解层级遇困

需要对比不同小波基与分解层级下图像去噪方法的性能，但现有代码仅支持单层小波分解，尝试使用wavedec2和waverec2实现多层级分解时，频繁遇到输入/输出参数数量错误。现有可运行的单层去噪代码如下：

% Load dataset
original = imread('brain_mri_transversal_t2_001.jpg'); % original image filename
noisy = imread('mri_image_noisy.jpg'); % Load noisy image

% Apply wavelet transform
[LL, LH, HL, HH] = dwt2(noisy, 'db4');

% Define the threshold value
t = 200;                  % Adjust the threshold value as needed
%t = rbe*sqrt(2*log(1000));      % tuni VisuShrink

% Define the function for adaptive thresholding
function coeff_thresh = adaptive_thresholding(coeff, threshold, rbe)
    if threshold == rbe
        % Adaptive Hard Thresholding for t = rbe
        coeff_thresh = zeros(size(coeff));
        coeff_thresh(coeff < -threshold) = coeff(coeff < -threshold);
        coeff_thresh(abs(coeff) <= threshold) = ...
            rbe*(exp((-coeff(abs(coeff) <= threshold).^2)/(2*(rbe^2))-1/2)) - rbe*(exp((-0)/(2*(rbe^2))-1/2));
        coeff_thresh(coeff > threshold) = coeff(coeff > threshold);
    elseif threshold > rbe
        % IMPROVED AGGD for t > rbe
        coeff_thresh = zeros(size(coeff));
        coeff_thresh(coeff < -threshold) = coeff(coeff < -threshold).*(1./(1+exp(coeff(coeff < -threshold)+threshold))-threshold/2);
        coeff_thresh(abs(coeff) <= threshold) = ...
            rbe*(exp((-coeff(abs(coeff) <= threshold).^2)/(2*(rbe^2))-1/2)) - rbe*(exp((-0)/(2*(rbe^2))-1/2));
        coeff_thresh(coeff > threshold) = coeff(coeff > threshold).*(1./(1+exp(coeff(coeff > threshold)+threshold))+threshold/2);
    end
end

% Robust median estimator 
rbe = median(abs(HH(:)))/0.6745;

% Threshold the wavelet coefficients
HH_thresh = adaptive_thresholding(HH, t, rbe);
LH_thresh = adaptive_thresholding(LH, t, rbe);
HL_thresh = adaptive_thresholding(HL, t, rbe);

% Inverse wavelet transform
reconstructed = idwt2(LL, LH_thresh, HL_thresh, HH_thresh, 'db4');

以下是重构后的代码，支持自定义分解层级和小波基，同时集成PSNR、SSIM性能指标计算，可直接用于不同参数组合的去噪效果对比：

% 加载图像
original = imread('brain_mri_transversal_t2_001.jpg');
noisy = imread('mri_image_noisy.jpg');
% 彩色图像转灰度图（若需处理彩色可扩展为多通道）
if size(original,3) == 3
    original = rgb2gray(original);
    noisy = rgb2gray(noisy);
end
original = double(original);
noisy = double(noisy);

% 定义待对比的参数组合
wavelet_list = {'db4', 'sym8', 'coif5', 'haar'}; % 常用小波基示例
level_list = [1, 2, 3, 4]; % 分解层级范围
threshold = 200; % 阈值参数

% 初始化性能结果表
results = table('VariableNames', {'Wavelet', 'Level', 'PSNR', 'SSIM'});

% 遍历所有参数组合进行去噪与评估
for wavelet_idx = 1:length(wavelet_list)
    current_wavelet = wavelet_list{wavelet_idx};
    for level_idx = 1:length(level_list)
        current_level = level_list(level_idx);
        
        % 多层级二维小波分解：C为系数向量，S为各层级尺寸矩阵
        [C, S] = wavedec2(noisy, current_level, current_wavelet);
        
        % 用最高层细节系数计算鲁棒中位数估计器rbe
        hh_coeff = detcoef2('hh', C, S, current_level);
        rbe = median(abs(hh_coeff(:)))/0.6745;
        
        % 逐层处理细节系数（LH、HL、HH）
        for l = 1:current_level
            % 提取当前层级的三个细节系数
            lh_coeff = detcoef2('lh', C, S, l);
            hl_coeff = detcoef2('hl', C, S, l);
            hh_coeff = detcoef2('hh', C, S, l);
            
            % 应用自定义自适应阈值处理
            lh_thresh = adaptive_thresholding(lh_coeff, threshold, rbe);
            hl_thresh = adaptive_thresholding(hl_coeff, threshold, rbe);
            hh_thresh = adaptive_thresholding(hh_coeff, threshold, rbe);
            
            % 将处理后的系数替换回C向量（注意wavedec2的系数存储顺序）
            C = wrev(C);
            C = [wrev(lh_thresh(:)); wrev(hl_thresh(:)); wrev(hh_thresh(:)); C(1+numel(lh_thresh)*3:end)];
            C = wrev(C);
        end
        
        % 重构去噪后的图像
        reconstructed = waverec2(C, S, current_wavelet);
        % 将像素值截断到图像合法范围[0,255]并转换为uint8格式
        reconstructed = uint8(clamp(reconstructed, 0, 255));
        
        % 计算性能指标：PSNR（峰值信噪比）、SSIM（结构相似性）
        psnr_val = psnr(reconstructed, uint8(original));
        ssim_val = ssim(reconstructed, uint8(original));
        
        % 存储当前参数组合的性能结果
        results = [results; {current_wavelet, current_level, psnr_val, ssim_val}];
    end
end

% 打印对比结果
disp('不同小波基与分解层级的去噪性能对比：');
disp(results);

% 自定义自适应阈值处理函数
function coeff_thresh = adaptive_thresholding(coeff, threshold, rbe)
    coeff_thresh = zeros(size(coeff));
    if threshold == rbe
        % 自适应硬阈值逻辑
        coeff_thresh(coeff < -threshold) = coeff(coeff < -threshold);
        mid_idx = abs(coeff) <= threshold;
        coeff_thresh(mid_idx) = rbe*(exp((-coeff(mid_idx).^2)/(2*(rbe^2)) - 1/2)) - rbe*(exp(-1/2));
        coeff_thresh(coeff > threshold) = coeff(coeff > threshold);
    elseif threshold > rbe
        % 改进AGGD阈值逻辑
        coeff_thresh(coeff < -threshold) = coeff(coeff < -threshold).*(1./(1+exp(coeff(coeff < -threshold)+threshold)) - threshold/2);
        mid_idx = abs(coeff) <= threshold;
        coeff_thresh(mid_idx) = rbe*(exp((-coeff(mid_idx).^2)/(2*(rbe^2)) - 1/2)) - rbe*(exp(-1/2));
        coeff_thresh(coeff > threshold) = coeff(coeff > threshold).*(1./(1+exp(coeff(coeff > threshold)+threshold)) + threshold/2);
    end
end

% 数值截断函数：确保像素值在合法范围内
function out = clamp(in, min_val, max_val)
    out = max(min(in, max_val), min_val);
end

关键说明 #

• 多层级分解处理：wavedec2返回的系数向量C按从粗到细的顺序存储，通过wrev反转后替换细节系数，再反转回去保证顺序正确，避免参数匹配错误。
• 参数扩展性：可直接修改wavelet_list和level_list添加更多待对比的小波基或分解层级。
• 图像格式兼容：自动处理彩色转灰度，重构后截断像素值到0-255范围，避免图像显示异常。
• 性能量化：使用MATLAB内置的psnr和ssim函数（若版本不支持，可自行实现对应公式），直观对比不同参数的去噪效果。

内容的提问来源于stack exchange，提问作者Mikasa Kirigiri