基于MATLAB的NSCT（非下采样轮廓波变换）实现

基于MATLAB的NSCT（非下采样轮廓波变换）实现代码及操作说明

---

一、环境配置与工具箱安装

1. 工具箱下载

   • 从CSDN文库或MATLAB File Exchange下载NSCT工具箱（如Nonsubsampled Contourlet Toolbox）。
   • 解压后添加到MATLAB路径：addpath('NSCT_Toolbox\')。

2. 关键函数说明

   • nsctdec: NSCT分解函数，支持多级多方向分解。
   • nsctrec: NSCT重构函数，根据分解系数重建图像。
   • wthcoefsp: 阈值处理函数，用于去噪。

---

二、NSCT分解与重构示例代码

% 读取图像并转换为灰度
img = imread('lena.png');
gray_img = rgb2gray(img);

% NSCT分解参数设置
levels = 3;       % 分解层数
directions = 8;   % 每层方向数（常用4/8/16）

% 执行NSCT分解
[pyr, dfb] = nsctdec(gray_img, levels, directions);

% 显示分解结果（低频与高频子带）
figure;
subplot(2,2,1); imshow(pyr{1}, []); title('低频分量 (LL)');
for i = 2:length(pyr)
    subplot(2,2,i); imshow(pyr{i}, []); 
    title(sprintf('高频分量 (Level %d)', i-1));
end

% NSCT重构
reconstructed_img = nsctrec(pyr, dfb);

% 显示原始与重构图像对比
figure;
subplot(1,2,1); imshow(gray_img); title('原始图像');
subplot(1,2,2); imshow(reconstructed_img, []); title('NSCT重构图像');

---

三、NSCT去噪完整流程

% 读取含噪图像
noisy_img = imread('noisy_lena.png');
noisy_gray = rgb2gray(noisy_img);

% NSCT分解（3层，8方向）
J = nsctdec(noisy_gray, 3, 8);

% 阈值处理（软阈值）
denoised_bands = cell(size(J));
for band = 1:numel(J)
    denoised_bands{band} = wthcoefsp('s', J{band}, 0.3, 'sln'); % 阈值参数可调
end

% NSCT重构
denoised_img = nsctrec(denoised_bands);

% 计算PSNR评估效果
psnr_value = psnr(denoised_img, noisy_gray);
disp(['PSNR: ', num2str(psnr_value), ' dB']);

% 可视化对比
figure;
subplot(1,3,1); imshow(noisy_gray); title('含噪图像');
subplot(1,3,2); imshow(denoised_img, []); title('去噪结果');
subplot(1,3,3); imshow(noisy_gray - im2double(denoised_img)); title('残差');

---

四、关键参数优化策略

1. 分解层数（levels）

   • 建议层数：3-4层（过多会导致计算量激增）。
   • 示例：levels =表示逐层分解。

2. 方向数（directions）

   • 低频子带（Level 1）：通常不分解。
   • 高频子带：方向数逐层增加（如“）。

3. 阈值选择

   • 全局阈值：适用于均匀噪声（如wthcoefsp('g', J, 0.2)）。
   • 自适应阈值：根据子带能量动态调整（需自定义函数）。

---

五、应用扩展（图像融合）

% 图像融合示例（以多聚焦图像为例）
img1 = imread('focus1.png'); img2 = imread('focus2.png');
gray1 = rgb2gray(img1); gray2 = rgb2gray(img2);

% NSCT分解
[J1, dfb1] = nsctdec(gray1, 3, 8);
[J2, dfb2] = nsctdec(gray2, 3, 8);

% 融合规则（低频取平均，高频取绝对值大者）
fused_pyr = cell(size(J1));
for i = 1:numel(J1)
    if i == 1
        fused_pyr{i} = (J1{i} + J2{i}) / 2; % 低频融合
    else
        fused_pyr{i} = max(abs(J1{i}), abs(J2{i})); % 高频融合
    end
end

% 重构融合图像
fused_img = nsctrec(fused_pyr, dfb1);
imshow(fused_img);

---

六、常见问题解决

1. 函数未找到错误

   • 检查工具箱路径是否添加：which nsctdec应返回函数路径。
   • 重新编译C文件（如mex atrouse.c）。

2. 分解效率低

   • 使用GPU加速（需Parallel Computing Toolbox）：

     pyr = nsctdec_gpu(gray_img, levels, directions);
     

3. 边界伪影

   • 添加周期延拓：在分解前执行gray_img = padarray(gray_img, [32,32], 'replicate')。

---

七、参考文献与工具

• 源码：NSCT轮廓波变换程序 www.youwenfan.com/contentald/52707.html
• 工具包：NSCT Toolboxes on File Exchange（推荐Contourlet Toolbox）。
• 论文参考：Donoho的轮廓波变换原始论文及NSCT改进文献。