简介:双树复小波变换(DT-CWT)是一种面向信号处理和图像分析的高级多尺度工具,尤其擅长图像融合中的细节保留与多方向边缘提取。这份资源将相关算法封装为轻量工具箱,适合需要实现或验证DT-CWT的研究人员、工程师及高年级学生。压缩包共45个文件,体积仅83KB,包含30个m函数脚本、11个mat测试数据、3个asv自动备份文件及1个txt说明;m文件覆盖正/逆变换、方向滤波器组、平移不变性测试和一维/二维演示脚本,mat文件则提供了Lenna等经典样本及不同阶数的Q-shift滤波器系数,方便直接调用。目前已有371人学习/下载。借助该工具箱,可在几分钟内完成双树复小波分解与重构,观察不同方向子带的差异;也可将其集成进图像融合流程,与传统小波、轮廓波等方法对比,支撑科研实验或课程设计。
1. 双树复小波变换工具箱:图像融合里的平移不变性与方向选择性
做医学图像融合时,最折磨人的往往不是融合规则,而是变换伪影。普通离散小波变换对平移敏感,配准误差稍大,融合边缘就出现振铃;方向选择性也只有三个方向,应对斜向纹理无能为力。DT-CWT用两棵并行小波树构造复数系数,同时获得平移不变性和六方向选择性,计算量只增加一倍,细节保留却明显更好。dtcwt_toolbox4_3是Kingsbury小组早期的MATLAB工具箱,内置一维/二维正反变换、Q-shift滤波器组、平移测试脚本与标准测试图,至今仍是论文复现的首选参考实现。想在图像融合里验证DT-CWT,它能让你跳过滤波器设计,直接进入算法实现。
2. DT-CWT 核心原理解读:两棵树换来平移不变性与六方向分解
2.1 双树结构与复数小波系数的由来
DT-CWT 之所以叫"双树",是因为它并行维护两棵独立的离散小波变换树:一棵称为实部树,一棵称为虚部树。两棵树使用相互关联但并不相同的滤波器组,重点在于虚部树相对实部树的采样位置错开半个采样间隔。这个设计带来的直接后果是,实部树的小波函数和虚部树的小波函数构成一对近似希尔伯特变换对,于是任意一维信号经过两棵树分解后,可以在每个尺度、每个位置得到一对系数,组合成复数小波系数:实部来自实部树,虚部来自虚部树。
平移不变性的来源就在这种互补采样结构里。普通DWT每一层只保留一个抽样相位,输入信号平移一个样本后,系数会在不同子带之间重新分配,重构结果跟着跳变。DT-CWT 的实部树以整数位置采样,虚部树以半整数位置采样,当输入平移一个样本时,两棵树系数变化方向相反,合成复数后模值几乎不变,只有相位发生线性旋转。这个相位旋转正好对应局部位移,也为后面的方向选择性提供了数学基础。
在 dtcwt_toolbox4_3 中,一维变换入口是 dtwavexfm.m,先并行调用两组分析滤波器,再用 q2c.m 把两棵树输出按奇偶抽样交叉组合成复数系数;二维变换 dtwavexfm2.m 在行方向和列方向各做一次同类操作,每一层输出一个低频近似子带和一组方向高频子带。逆变换 dtwaveifm.m 与 dtwaveifm2.m 把复数系数拆回实部、虚部,分别送入两棵树的综合滤波器再求和。整个过程要求分析滤波器组与综合滤波器组满足完全重构条件,否则重建图像会出现系统性的灰阶失真。
与普通DWT相比,DT-CWT的计算量大约只多一倍:两棵树并行,每棵树都是标准滤波器组;二维情况下,行和列各做两次滤波,总复杂度仍然是O(N²)量级,不会因为复数表示而变成无法接受的高阶代价。这也是它能在图像融合场景下替代DWT的重要原因——性能收益明显,成本可控。
2.2 方向选择性:为什么是 6 个方向而不是 3 个
二维 DWT 沿行、列各做一次一维分解,产生 LL、LH、HL、HH 四个区域,方向信息只能区分水平、垂直和对角三个方向。DT-CWT 因为每层都有实部和虚部两套小波,在二维情况下行方向和列方向分别选实部或虚部,可以得到多种组合;实部与虚部之间的希尔伯特关系让真正独立的方向信息收敛为 6 个:约 ±15°、±45°、±75°。这组角度来自 Kingsbury 对滤波器组相位响应的推导,在工具箱的二维分解结果中,高频元胞数组的第三维索引依次对应:
| 方向子带索引 | 对应角度 | 主要捕捉的结构 |
|---|---|---|
| 1 | +15° | 接近水平的斜边、细长纹理 |
| 2 | -15° | 接近水平的反向斜边 |
| 3 | +45° | 主对角线边缘 |
| 4 | -45° | 副对角线边缘 |
| 5 | +75° | 接近垂直的斜边 |
| 6 | -75° | 接近垂直的反向斜边 |
与 DWT 相比,DT-CWT 的关键改进在于水平、垂直和斜向边缘不再混在同一个子带里,而是被拆到不同方向通道。融合时可以针对某个方向单独设置权重,这在医学图像融合里特别实用,因为骨骼边缘与软组织纹理往往落在不同方向子带上。
2.3 从 dtwavexfm2.m 看分解输出的数据结构
实际使用这个工具箱,第一件要弄清的事是输出格式。二维正变换的标准调用方式是:
[Yl, Yh] = dtwavexfm2(X, nlevels, 'near_sym_a', 'qshift_a');X 为输入灰度图像,要求尺寸为偶数,内部默认按反射边界扩展。Yl 是低频近似系数矩阵,大小约为输入图像在最后一层缩放后的尺寸;Yh 是长度为 nlevels 的元胞数组,第 level 层是一个 H×W×6 的复数数组,第三维依次对应当前层的 6 个方向子带。
参数说明:nlevels 是分解层数,一般取 3 到 5。层数越多,低频子带越粗糙,方向信息越抽象,边界效应也会累积。'near_sym_a' 是第一层实部树和虚部树共用的短滤波器组,'qshift_a' 是后续层使用的 Q-shift 滤波器组。第一层不用 Q-shift 是因为 Q-shift 滤波器要求信号预先按半采样间隔对齐,直接处理原始像素会放大边界偏差;near_sym 这类接近对称的短滤波器更适应原始像素网格。Q-shift 名字里的 Q 代表 quarter-sample,即四分之一采样延迟,专门用来把两棵树之间的相位差稳定在 90°,保证实部虚部构成解析对。
逆变换的输入输出关系完全对称:
recon = dtwaveifm2(Yl, Yh, 'near_sym_a', 'qshift_a');如果分解和重构用的滤波器组参数不一致,输出会出现明显棋盘格纹理。如果输入不是 double 类型,uint8 图像在卷积时会被截断到 255,重构误差会放大到肉眼可见的程度。所以进入变换前先转 double,是使用这个工具箱的第一个习惯动作。
3. 动手复现:部署 dtcwt_toolbox4_3 并跑通一维/二维变换实验
3.1 工具箱目录结构与关键文件定位
解压 dtcwt_toolbox4_3.rar 后,目录里没有复杂的工程结构,而是平铺的 .m 源文件与 .mat 数据文件。按功能可以分成几组:
| 类别 | 文件 | 作用 |
|---|---|---|
| 一维变换 | dtwavexfm.m / dtwaveifm.m | 一维 DT-CWT 分解与重构 |
| 二维变换 | dtwavexfm2.m / dtwaveifm2.m | 二维分解与重构,图像融合主入口 |
| 滤波器操作 | colfilter.m / colifilt.m / coldfilt.m / coliwtfilt.m | 列方向滤波,带 i 的为逆或隔点采样版本 |
| 系数组合 | q2c.m | 把多路滤波输出重组为复数小波系数 |
| 边界处理 | reflect.m | 反射式边界扩展,减少边缘效应 |
| 测试脚本 | shift_test_1D.m / shift_test_2D.m / shift_test_2DAA.m / shiftmovie.m | 验证平移不变性的实验程序 |
| 测试数据 | lenna.mat / qshift_a~d.mat / near_sym_a.mat / near_sym_b.mat / antonini.mat / legall.mat | 标准测试图与滤波器组系数 |
| 可视化 | cimage5.m / setfig.m / SETTITLE.M | 复数系数彩色显示与图像设置 |
建议把整个目录加入 MATLAB 路径,而不是 cd 到目录里运行。工具函数之间以相对名字互相调用,加入路径后在任何工作目录都能直接调用 dtwavexfm2。
注意:只添加 dtcwt_toolbox4_3 解压后的文件夹本身,不要连同外层压缩包目录一起添加,避免与 MATLAB 自带同名函数冲突。
3.2 一维平移不变性对照:shift_test_1D.m 到底在测什么
直接运行 shift_test_1D.m 会生成一个对比图,左列是普通 DWT 在不同平移量下的重构结果,右列是 DT-CWT 的结果,每一行对应信号平移 0、1、2、3 个样本。观察重点不是重构曲线本身,而是同一列里不同平移量之间的差异。DWT 一侧会因为输入平移导致小波系数在子带间重新分配,重构曲线出现明显跳变;DT-CWT 一侧由于两棵树互补采样,不同平移量下的重构曲线几乎重叠。
如果不想跑完整脚本,可以手动复现这个实验。下面用工具箱自带函数完成一维测试:
x = zeros(256, 1); x(64:68) = 1; % 构造窄脉冲,接近理想冲击信号 for shift = 0:3 xs = circshift(x, shift); % 循环平移 0~3 个样本 [a, d] = dtwavexfm(xs, 3, 'near_sym_a', 'qshift_a'); xr = dtwaveifm(a, d, 'near_sym_a', 'qshift_a'); err(shift + 1) = max(abs(xr - xs)); end disp(err);逻辑说明:dtwavexfm 对平移后的信号做三层分解,低频系数 a 和高频元胞 d 一起送进 dtwaveifm 重构,然后计算重构信号与原始信号的最大绝对误差。如果 err 四个分量几乎相同,说明平移没有显著改变重构结果,即平移不变性成立。作为对照,把 dtwavexfm/dtwaveifm 换成工具箱里 wavexfm/waveifm 的对应调用,再跑一遍同样的循环,DWT 侧的误差波动通常会比 DT-CWT 大一到两个数量级。
参数说明:分解层数取 3 对一维信号足够。circshift 是循环平移,避免了边界补零带来的额外误差,比手动索引平移更稳。
工具箱里还有一个 shift_test_2DAA.m,区别在于它使用 Allpass 滤波器链做分数平移(包括 0.25、0.5 像素),测试的是亚像素平移敏感性。这个脚本运行较慢,但能给出更接近真实配准残差的评估,融合前值得跑一次。
3.3 二维图像分解:从 lenna.mat 到 dtwavexfm2/dtwaveifm2
二维场景我习惯先用工具箱自带的 lenna.mat 做完整重构验证,环境没问题再换成自己的数据。
load lenna.mat; % 载入工具箱自带测试图 im = double(lenna); % uint8 转 double,避免卷积截断 nlevels = 4; [Yl, Yh] = dtwavexfm2(im, nlevels, 'near_sym_a', 'qshift_a'); recon = dtwaveifm2(Yl, Yh, 'near_sym_a', 'qshift_a'); err = max(abs(recon(:) - im(:))); fprintf('最大重构误差: %.4g\n', err);逻辑说明:dtwavexfm2 把 lenna 分解为 4 层,Yl 是低频近似,Yh 是 4×6 个方向子带,然后 dtwaveifm2 用同样参数重构。err 在 1e-10 量级时说明分析滤波器与综合滤波器配对正确,变换可逆,后续融合时对系数的修改才能稳定映射回图像域。
这里有三个容易踩的坑。第一,lenna.mat 里保存的是 0~255 的 uint8 图像,不转 double 的话卷积结果被截断,重构误差会从 1e-10 跳到 0.1 量级。第二,nlevels 不能超过 log2(min(size(im))),否则最内层子带尺寸小于滤波器长度,直接报维度不匹配。第三,重构时滤波器组必须与分解时完全一致,混用 near_sym_a 和 near_sym_b 会让图像出现棋盘格状高频噪声。
注意:如果 load 之后工作区没有 lenna 变量,用 whos 查看实际变量名,不同版本保存的变量名可能不同。
4. 图像融合实战:基于 DT-CWT 的高低频合并规则与参数调优
4.1 融合为什么选 DT-CWT 而不是 DWT
图像融合的目标是把同一场景下不同传感器或不同成像条件的多张图像合并为一张信息更完整的图。CT 对骨骼等硬组织清晰,MRI 对软组织清晰,融合目标是让两者同时可见,这是医学图像融合的典型诉求。DWT 在融合里的主要问题是平移敏感和方向性不足:平移敏感导致源图像配准残差被放大成融合伪影,方向性不足则让斜向边缘在合并时被平均操作削弱。
DT-CWT 的近似平移不变性让融合规则对配准误差更鲁棒,源图像即使存在亚像素级错位,也不会在融合结果里产生明显振铃带。六方向分解则让不同倾斜角度的边缘落在各自对应子带,合并时可以分别比较、单独选择,而不是把斜边信息混在水平垂直子带里做一刀切决策。这让 DT-CWT 在图像融合里长期作为经典对比基准,也是这个工具箱存在的主要意义。
另一个典型场景是多聚焦图像融合:同一场景下不同焦距拍两张图,一张前景清晰、一张背景清晰,融合目标是得到全景清晰图像。这类图像不存在传感器差异,但深度不连续导致的边缘位移更明显,DT-CWT 的六方向分解能分别处理前景轮廓和背景纹理,避免 DWT 中常见的振铃叠加。
4.2 系数合并策略:低频平均、高频模极大
融合框架固定为:源图像 → DT-CWT 分解 → 系数合并 → 逆变换 → 融合图像。合并规则决定融合质量,工具箱本身只提供变换能力,合并代码需要自己写。最常见的策略是低频取平均、高频取模极大值:
function F = dtcwt_fuse(A, B, nlevels) % A、B: 已配准且尺寸相同的灰度图像 % nlevels: 分解层数,医学图像推荐 4~5 if ~isa(A, 'double'), A = double(A); end if ~isa(B, 'double'), B = double(B); end [Al, Ah] = dtwavexfm2(A, nlevels, 'near_sym_a', 'qshift_a'); [Bl, Bh] = dtwavexfm2(B, nlevels, 'near_sym_a', 'qshift_a'); % 低频系数取平均,保留两幅图共有的亮度结构 Fl = 0.5 * (Al + Bl); % 高频系数逐层逐方向取模值大的 for level = 1:nlevels for dir = 1:6 mA = abs(Ah{level}(:, :, dir)); mB = abs(Bh{level}(:, :, dir)); mask = mA >= mB; Fh{level}(:, :, dir) = Ah{level}(:, :, dir) .* mask + ... Bh{level}(:, :, dir) .* (1 - mask); end end F = dtwaveifm2(Fl, Fh, 'near_sym_a', 'qshift_a'); F = uint8(F); % 转回显示范围 end逻辑说明:低频系数表示图像的大尺度亮度分布,CT 和 MRI 的低频信息来源不同但都比较平滑,直接平均可以在保留公共解剖结构的同时抑制亮度偏移。高频系数是复数,实部和虚部都携带纹理信息,直接比较实部没有意义,必须用 abs() 求模值。模值大表示当前位置能量更强,边缘或纹理更显著,所以保留模值较大的系数。mask 是逻辑矩阵,参与复数乘法时自动按 0/1 数值处理。
参数说明:nlevels 取 3 时计算快、方向信息粗糙,取 5 时细节更好,但底层子带尺寸只有原始图的 1/32,靠近边界的位置会积累较多边界效应。对 512×512 的医学图像,4 层是经验折中。两幅源图尺寸必须完全一致,否则 Yh 元胞维度对不上,合并阶段直接报错。F 最后转 uint8 的前提是输入 A、B 都在 0~255 范围,如果输入是 0~1 的归一化数据,去掉这行即可。
如果需要更精细的控制,可以把高频合并从二值选择改成加权组合,权重由两幅图系数的模值比例决定:w1 = mA ./ (mA + mB + eps),w2 = 1 - w1,然后 Fh = w1 .* Ah + w2 .* Bh。这个做法在模值接近的区域产生平滑过渡,融合结果比硬阈值更自然,医学图像融合里经常用。
4.3 融合效果验证与参数调优建议
融合结果不能只凭肉眼判断,我一般同时看三个指标。重构一致性用结构相似度 SSIM 衡量融合图与各源图的结构相似程度;边缘保持度用方向信息保留率 QAB/F 量化边缘和纹理的保留比例;灰度分布则看直方图是否有双峰断层或饱和。MATLAB 的图像处理工具箱直接提供 ssim:
ssim_A = ssim(F, A); % 靠近1说明与源图A结构一致 ssim_B = ssim(F, B); fprintf('SSIM_A = %.4f, SSIM_B = %.4f\n', ssim_A, ssim_B);逻辑说明:这里传入的 F、A、B 都应是 uint8 或 0~255 的 double,ssim 函数内部会统一处理动态范围。SSIM_A 和 SSIM_B 分别反映融合图保留两幅源图结构信息的能力,两者都高说明融合没有明显偏袒某一幅输入。
调参时先固定滤波器组为 near_sym_a + qshift_a,只改 nlevels 对比 3/4/5 层结果;然后固定层数,把滤波器组换成 near_sym_b 或 qshift_b,观察斜向纹理变化。低频平均策略在 CT/MRI 融合里比较稳定,但在多聚焦图像上会让纹理区域变糊,需要改成低频模极大或加权平均。
医学图像融合里还有一个实用技巧:先对两幅源图做直方图匹配,让亮度范围对齐,再做 DT-CWT 融合。否则低频平均会在骨骼与软组织交界处产生一条类似水肿的过渡带,实际是亮度失配造成的伪影。
5. 换滤波器组与看系数:把工具箱用到非标准场景的几个技巧
5.1 用 qshift_b / near_sym_b 调整滤波器响应
工具箱提供了 qshift_a 到 qshift_d、near_sym_a、near_sym_b 等滤波器组文件。qshift 系列都是 Q-shift 设计,区别在长度和频率选择性:qshift_a 长度较短、边界效应小,适合中小尺寸图像;qshift_d 更长、频率选择性更好,但相位线性稍差。替换方法很简单,把 dtwavexfm2 和 dtwaveifm2 的第三、第四参数改掉即可。判断时机:图像纹理呈周期性,比如织物、栅栏或医学影像里的骨小梁结构,用更长的 qshift_d 更容易把频率邻近的纹理分开;图像尺寸小且边缘锐利,比如 128×128 的局部病灶区域,用 near_sym_a 可以减少跨边界扩散。
5.2 借助 shift_test_2D.m 排查配准误差影响
正式融合前,先用 shift_test_2D.m 摸清当前参数下的平移敏感度。在命令窗口直接输入 shift_test_2D 即可运行。脚本对同一幅二维图像做多个整数像素平移,分别做分解再重构,在图上标出最大误差分布。如果差异图出现明显条状伪影,说明当前滤波器组在该尺寸下边界效应过大,应换更短滤波器或降低分解层数。实际任务中我会先把源图像相对平移控制在整数像素内,再用这个脚本确认重构误差在可接受范围,才开始正式融合,避免把配准残差造成的伪影误判成融合规则的问题。
5.3 用 cimage5.m 看复数系数的振幅与相位
高频子带系数是复数,直接 imshow 只能显示实部,丢掉相位信息。cimage5.m 是工具箱自带的复数图像显示函数,把振幅映射为亮度、相位映射为色相,一张图里同时看能量分布和相位结构:
figure; cimage5(Yh{2}(:, :, 3)); % 查看第2层、第3个方向子带 colorbar;相位信息对边缘定位很有用:同一方向子带里相位突变的位置往往对应亚像素级边缘,振幅峰值则对应强边缘中心。处理医学图像融合时,我习惯先对两幅源图分别显示几个关键方向子带的 cimage5,确认它们在哪些方向子带上差异最大,再决定是否需要为特定方向增加融合权重。这种方向层面的精细控制,是 DWT 融合里很难做到的,也是把 DT-CWT 工具箱从简单对比实验推向实际项目时最值得花时间的一步。
本文还有配套的精品资源,点击获取