conturlet变换-nsct_toolbox⼯具箱源码简要分析
因为需要编写关于contourlet的相关代码，苦于⽹上没有太多的相关资料，⽆奈，硬着头⽪强⾏啃⾷nsct_toolbox的源码，现总结出⼀些经验和⾃⼰的理解，供苦×的⼴⼤学⼦参考，⾥⾯还有很多疑问点，希望看到该篇⽂章的朋友不吝赐教，有错误之处还请指正。
y = nsctdec(x, levels, [dfilt, pfilt] )
这就是⾮下采样（也可选择采样滤波器）的contourlet分解的主函数
y = nsctdec(x, levels, [dfilt, pfilt] )
x：图像的输⼊矩阵，双精度数据类型
levels:⽅向滤波器组分解层数向量，注意是⼀组向量值
dfilt:⽅向滤波器组
pfilt：塔式分解的滤波器组
out：输出的结果是⼀组向量集合，集合第⼀个向量是低频分量的⼆维数组，接下来是从⼩到⼤排列⾼频⼦带。
⾸先是⼀些滤波器的准备：
filters = cell(4) ;//初始化⼀个数组存放滤波器
[h1, h2] = dfilters(dfilt, 'd');//得到⽅向⼆维滤波器组，⼀低⼀⾼
//得到⽐例尺度
h1 = h1./sqrt(2) ;
h2 = h2./sqrt(2) ;
//这个函数⽤来进⼀步⽣成滤波器，参数c表⽰column列⽅向
filters{1} = modulate2(h1, 'c');
filters{2} = modulate2(h2, 'c');
//通过线性组合，由原来⼆⽅向调制四⽅向
[filters{3}, filters{4}] = parafilters( h1, h2 ) ;
//⽣成塔式分解的滤波器
[h1, h2, g1, g2] = atrousfilters(pfilt);
滤波器准备完毕，主要是塔式分解来进⾏层数切割，⽅向滤波器对每⼀层的⾼频进⾏处理下⾯来看⼀下索引控制
//获得向量长度表⽰层数
clevels = length( levels ) ;
//这个是表⽰有多少层，+1是为了存放多出来的低频
nIndex = clevels + 1 ;
% 初始化输出向量数组
y = cell(1, nIndex) ;
关键代码如下：
%⾦字塔分解
for i= 1 : clevels
%该函数稍后介绍，是计算⾦字塔分解的函数，xlo是较粗尺度（下⼀轮的迭代分解），xhi是⾼频    [xlo, xhi] = nsfbdec(x, h1, h2, i-1) ;
if levels(nIndex-1) > 0
%对⾼频进⾏⽅向分解，这⾥⼜冒出⼀个进⼀步分解函数，第三个参数是剩余分解的层数
xhi_dir = nsdfbdec(xhi, filters, levels(nIndex-1));
y{nIndex}=xhi_dir ;%创建该层的⼦带向量
else
% ⼩于0代表分解结束，直接输出低频结果
y{nIndex}=xhi ;
end
% 更新层数索引
nIndex = nIndex - 1 ;
%进⼊新的循环啦：）
x = xlo ;
end
function [h0, h1] = dfilters(fname, type)
⽣成⽅向⼆维滤波器组，
fnmae：可选的滤波器名字，具体在代码⾥有注释
type:d或者c表⽰分解或者重构
代码就不关⼼了，数学问题了。只要知道这个函数提供了低通和⾼通滤波器就⾏。
同样的还有⼏个滤波器构造函数就不深⼊探索了。
function [y0, y1] = nsfbdec( x, h0, h1, lev )
这是⾮下采样的的⾦字塔变换。
x：需要被分解的图像，第⼀次就是原图，第⼆次到最后⼀次都是上⼀次的低频分量。
h0,h1:
function [y0, y1] = nsfbdec( x, h0, h1, lev )
if lev ~= 0
I2 = eye(2);
shift = -2^(lev-1)*[1,1] + 2; L=2^lev;
y0 = atrousc(symext(x,upsample2df(h0,lev),shift),h0,I2 * L);
y1 = atrousc(symext(x,upsample2df(h1,lev),shift),h1,I2 * L);
else
//到了第⼀层直接上⼩波变换（疑问：为什么contourlet第⼀层不能多⼏个⽅向呢，要⽤⼩波四⽅向？）
shift = [1, 1]; % delay compensation
//⼩波函数，注意⼯具箱⾥没有这个函数，⼤佬是⽤c写的，要在控制台编译相关c mex filename
y0 =  conv2(symext(x,h0,shift),h0,'valid');
y1 =  conv2(symext(x,h1,shift),h1,'valid');
end
function y = nsdfbdec( x, dfilter, clevels )
这是第⼆个主要函数，contourlet的核⼼，对⾼频进⾏⽅向分解了。先拿个⽼图镇贴。
在图上我们看到分解是逐级变多的。
X：输⼊图像，这⾥当然就是塔式分解下来的⾼频
dfilter:⼀个字符串，⽅向过滤器名称。矩阵单元，包括两个⽅向滤波器和⼋个平⾏四边形过滤器。刚刚准备的过滤器派上⽤场了。
⾸先得细说⼀下传⼊进来的索引控制，当初关于层数分解是由图像⼤⼩或者需求决定的，级数分解在传⼊进来的数组⾥都是同⼀个数，下⾯要分类讨论级数分解
这部分源码稍微有点复杂了
Input check和滤波器准备不说了
q1 = [1, -1; 1, 1];//先定义了⼀个采样矩阵
//分解级数为1， NSSFBDEC⽤⼀个双通道⾮⼦采样滤波器组对图像X进⾏卷积，⼀共就两个分⽀的输出modulate
if clevels == 1
[y{1}, y{2}] = nssfbdec( x, k1, k2 ) ;
//分解级数为2
//先把⼀级分解了
[x1, x2] = nssfbdec( x, k1, k2 ) ;
//根据⼀级分解下来的再次分解，所以级数索引当初设定是要取2的对数
[y{1}, y{2}] = nssfbdec( x1, k1, k2, q1 ) ;
[y{3}, y{4}] = nssfbdec( x2, k1, k2, q1 ) ;
//三级往上⼤概也是这么个意思。
到这⾥也基本上稀⾥糊涂的算是把主要函数把握了。等等，有啥⽤啊，分解的这么开⼼，都分解的稀巴烂了，⼜是纵向分解，⼜是横向分解，⼲嘛搞这么累。我们其实分解下来的⽬的就是在不同尺度不同⽅向下对频域进⾏处理，来达到我们的效果，⾄于怎么处理那就是要实验研究的地⽅了，处理完毕，直接交给我们的重构函数，⼀键⽣成！！让开让开，重构函数来了。
function x = nsscrec(y, dfilt, pfilt)
y：就是分解下了的结果，处理的时候注意不破坏原有的数据结构就⾏
还有两个参数同原来的分解，怎么分解怎么还原。
输出的就是⼀个图像的矩阵。
for i=1:n
% Process the detail subbands
if iscell( y{i+1} )
//各个⽅向的重建
xhi = nsdfbrec( y{i+1}, filters );
else
//如果没有⽅向分解，直接原路返回
xhi = y{i+1};
end
%⾦字塔重建
x = nsfbrec(xlo, xhi, g1, g2, nIndex);
%为下⼀次重建准备
xlo = x ;//这⼀层重建完成的就是下次的低频，具体看上⾯的图就理解了
nIndex = nIndex -1;
end

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。