下载文档原格式
基于双树复小波变换的图像去噪的开题报告一、选题背景随着数字图像技术的发展，图像噪声问题成为了图像处理领域中的一个重要问题，尤其是在数字信号传输、数字摄影、数字视频压缩等领域中，图像噪声问题更是不可避免。
因此，图像去噪技术成为了图像处理领域中的研究热点问题之一。
在现代图像处理技术中，基于小波变换的去噪技术已经被广泛应用。
小波变换方法可以将图像数据通过一定的方法进行变换，使得信号特征得到更好的表现，同时还可以减小信号噪声的影响。
但是在实际应用中，小波变换的去噪方法也存在一定的缺陷，如对于高斯噪声和脉冲噪声效果并不理想等。
为了解决小波变换方法的缺陷，研究人员提出了双树复小波变换去噪方法。
双树复小波变换是小波变换的一种扩展方法，其能够提供更好的信号描述力和更准确的信号特征表示。
因此，双树复小波变换在图像处理领域中具有较高的研究价值和应用前景。
二、研究目的本研究的目的是探究基于双树复小波变换的图像去噪方法，研究双树复小波变换的原理及实现方法，对其进行分析、优化和改进，提出一种基于双树复小波变换的图像去噪算法，并进行实验验证，以提高图像噪声去除的效率和精度。
三、研究内容1. 双树复小波变换的原理与实现方法的研究双树复小波变换是对小波变换的扩展方法，该方法采用双树结构和复数扩展，能够更好地描述信号特征。
本研究将学习双树复小波变换的原理，介绍其在信号处理中的优劣势，并实现其算法。
2. 基于双树复小波变换的图像去噪算法的设计在研究双树复小波变换的基础上，本研究将探讨其在图像处理中的应用，并设计一种基于双树复小波变换的图像去噪算法，该算法能够有效消除噪声影响并保持图像细节特征。
3. 基于算法的实验验证在算法设计完成后，本研究将进行实验验证。
该实验将采用一系列有噪图像，分别对比本算法与其他常见的基于小波变换的图像去噪算法，包括小波软阈值去噪算法、小波硬阈值去噪算法等，分析其去噪效果、处理速度等指标，并比较其优缺点。
四、研究意义1. 提高图像去噪的效率和精度，满足实际应用需要双树复小波变换能够提供更好的信号特征描述能力，因此在图像去噪中具有一定优势。
E HRStffI i N r A a«VT创新助手报告——主题分析报告创新助手平台提供北京万方软件股份有限公司Array2014-07-10报告目录报告核心要素 ....................................................... I.一、主题简介 ...................................................... 1.二、主题相关科研产出总体分析 (1)2.1文献总体产出统计 (1)2.2学术关注趋势分析 (2)三、主题相关科技论文产出分析 (2)3.1中文期刊论文................................................ 2.3.1.1近十年中文期刊论文分布列表 (2)3.1.2中文期刊论文增长趋势 (3)3.1.3发文较多期刊 (4)3.1.4发文较多的机构 (4)3.1.5发文较多的人物 (4)3.1.6核心期刊分布数量对比 (4)3.1.7最近相关中文期刊论文 (5)3.1.8被引较多的相关期刊论文 (6)3.2学位论文.................................................... 7.3.2.1近十年学位论文年代分布列表 (7)3.2.2学位论文增长趋势 (8)3.2.3硕博学位论文数量对比 (8)3.2.4发文较多的机构 (9)3.2.5发文较多的人物 (9)3.2.6最近相关学位论文 (9)3.3中文会议论文 (10)3.3.1近十年中文会议论文年代分布列表 (10)3.3.2中文会议论文增长趋势................................. 1.13.3.3中文会议论文主办单位分布............................. 1 13.3.4发文较多的机构....................................... 1.23.3.5发文较多的人物....................................... 1.23.3.6最近相关中文会议论文................................. 1.33.4外文期刊论文 (13)3.4.1近十年外文期刊论文年代分布列表 (13)3.4.2外文期刊论文增长趋势................................. 1.43.4.3最近相关外文期刊论文 (14)3.5外文会议论文 (15)。
基于双树复小波变换的图像检索新算法舒彬【摘要】本文研究了一种双树复小波变换的图像形状特征检索算法。
第一步，利用双树复小波对原图像进行变换，取各层的2个低频子图像，分别求出它们的模极大值的平均值和原图像的模极大值的平均值；第二步，将它们分别作为下一层的阈值，找出每层的边缘图像；第三步，通过计算不变矩之间的欧式距离，查询出与原图相似的图像。
实验结果表明，此算法的检索效率比canny算子的算法和Mallt小波变换的算法都高。
【期刊名称】《产业与科技论坛》【年(卷),期】2016(015)009【总页数】2页(P57-58)【关键词】双树复小波变换;模极大值;动态阈值;不变矩;图像检索【作者】舒彬【作者单位】陕西学前师范学院【正文语种】中文形状是图像最基本的特征之一。
提取形状特征对图像进行检索，是一种很重要的方法。
因为具有很强的多尺度分析能力[1]和非常良好的时频局域化性质，从而使双树复小波变换很适合检测目标图像的概貌和微小的细节。
表述形状特征的时候[2]，表示它们的描述符应该具备如下优点：旋转、平移和尺度不变形等。
尽管对于描述物体形状的方法有很多，可它们大多给的结论不详细，对边缘的表述不好理解。
基于此，本文研究了利用双树复小波变换对图像的形状特征进行检索的算法。
在信号分析中，如果用离散小波变换分析形状特征的话，非常微小的平移，便会引起小波系数在不同尺度上的很大差异，尽管形状相同的两个物体，产生的形状特征向量却不相同。
它缺乏方向选择性，弱化了其它方向的信息，从而丢失了一些重要的信息[3]。
为了克服这些问题，Kingsbury等人提出了双树复小波变换(Dual-tree Complex Wavelet Transform,DT-CWT)[4]。
双树复小波变换采用两个实小波可称为树A 和树B。
其中上部树A表示变换的实部，下部树B表示变换的虚部。
树A及树B 变换采用两个不同的滤波器集合，它们都符合重建条件。
假设h0(n)表示树A对应的滤波器组的低通滤波器对，h1(n)表示高通滤波器对，g0(n)表示树B对应的滤波器组的低通滤波器对，g1(n)表示高通滤波器对。
双树复小波变换
双树复小波变换（Dual-Tree Complex Wavelet Transforms, DT-CWT）是一种
为定位有用的（directional）信号分析中的关键技术，它利用相同的波形分解加
以微调，除了满足基于小波的特征分析需求外，还采取双树的结构，以确保理论预期的增益。
DT-CWT具有一种较强的完整性（coherence），可以有效地捕获复杂的
频率特征，即使是在低频域存在噪声的情况下仍可突出高频特征，因此受到了广泛的关注。
DT-CWT具有若干优势，如它有良好的分辨率，可以提取频带范围较窄的信号；具有几乎抗磁性的鲁棒性，能够具备良好的微弱信号抽取性能；另外它还具有优越的均衡性，可以很好地处理数据不均匀性，实现对复杂应用场景的有效建模；其次，它能够有效识别一时间周期分量，从而进一步提高定位能力。
作为一种基于实值小波变换的先进技术，DT-CWT非常适用于图像处理过程，
能够满足不同信号处理任务的要求，如图像压缩、图像分类、图像定位及图像剪裁等，为图像处理中高精度的实现奠定了基础。
DT-CWT在其他领域的应用也同样引
起了广泛的关注，如信号处理中定位低频信号、改善医学影像质量及构建语音自动识别系统等应用都能体现出其决定性的作用。
可以看出，DT-CWT不仅是一种简单有效的信号分析工具，而且对复杂信号处
理有着广泛的应用前景，未来估计可期预料其研究范围和领域的进一步扩大。
基于双密度双树复小波变换的目标检测和跟踪算法研究的开题报告一、选题背景目标检测和跟踪是计算机视觉领域的重要研究方向之一。
在实际应用中，它们广泛应用于智能监控、自动驾驶、无人机、机器人等领域。
目标检测主要是通过对图像或视频进行处理，找出其中的感兴趣目标，而目标跟踪则是跟踪目标在连续帧中的位置和运动轨迹。
传统的目标检测和跟踪方法通常采用特征提取和分类器技术，这些方法已经取得了不错的效果，但是在面对复杂场景、目标形变和遮挡等问题时存在一定的局限性。
为了克服这些问题并提高检测和跟踪的效果，近年来涌现出了一些基于深度学习、卷积神经网络等技术的新方法。
本课题将针对目标检测和跟踪中存在的问题，探讨一种基于双密度双树复小波变换的新方法。
该方法结合了小波变换和复小波变换的优势，可以提高图像处理的精度和速度，同时保持较好的图像局部特征。
二、研究内容本课题的主要研究内容包括以下方面：1.研究双密度双树复小波变换的理论和算法，了解其优势和局限性。
2.基于双密度双树复小波变换构建目标检测和跟踪模型，实现对感兴趣的目标的提取和跟踪。
3.对模型进行优化和测试，对比传统的目标检测和跟踪方法的效果，并对模型的性能进行评估和改进。
4.设计相关实验，验证方法的实用性和效果，并对实验结果进行分析和总结。
三、研究意义本课题采用基于双密度双树复小波变换的目标检测和跟踪方法，有以下研究意义：1.充分发掘小波变换和复小波变换的优势和特点，提高图像处理的精度和速度。
2.优化传统的目标检测和跟踪方法，提高检测和跟踪的成功率和鲁棒性。
3.扩展目标检测和跟踪的应用领域，使其在智能监控、自动驾驶、无人机、机器人等领域更具实用性。
四、研究方法本课题主要采用实验法和分析法进行研究，具体方法如下：1.收集本领域的相关文献和数据集，对相关技术进行深入了解和研究，为构建目标检测和跟踪模型提供理论支持。
2.设计双密度双树复小波变换的目标检测和跟踪模型，并优化模型框架和参数，提高检测和跟踪的准确性和鲁棒性。
---文档均为 word 文档，下载后可直接编辑使用亦可打印-要噪声抑制是任何图像处理任务的组成部分，噪声会显着降低图像质量，因此使观察者难以区分图像的细节，特别是在诊断检查中。
经过几十年的研究，已经提出了大量关于图像去噪的方法。
通过使用空间滤波或变换域滤波，可以减少图像中噪声的影响。
在变换域小波方法中，提供更好的去噪效果，同时保留像边缘那样的图像细节。
离散小波变换具有一些缺点，即由于缺乏移位不变性和较差的方向选择性，导致其在图像处理中的应用尚未确定。
为了克服这些缺点，使用了双树复数小波变换，其在传统的小波变换上提供了完美的重构。
它使用 2 个离散实小波变换；第一个离散实小波变换给出了变换的实部，而第二个离散实小波变换给出了变换的虚部。
它在二维和更高维度上有限的冗余几乎是不变和定向选择性的。
双树复数小波变换在图像去噪和增强等应用方面优于离散小波变换。
双树复数小波变换的优点之一是它可用于实现比二维离散小波变换方向更具选择性的二维小波变换。
二维双树复数小波在每个尺度上产生十二个子带，每一个都以不同的角度精确定位。
关键词：图像；去噪；双树复数小波；阈值ABSTRACTNoise suppression is an integral part of any image processing task. Noisesignificantly degrades the image quality and hence makes it difficult for the observer todiscriminate fine detail of the images especially in diagnostic examinations. Throughdecades of research, a lot of methods on image denoising have been proposed .Theeffect of noise in the images can be reduced by using either spatial filtering or transformdomain filtering. In transform domain, the wavelet method provides better denoisingeffect while preserving the details of images like edges. The Discrete WaveletTransform (DWT) has some disadvantages that undetermined its application in imageprocessing as lack of shift invariance and poor directional selectivity. In order toovercome these disadvantages Dual Tree Complex Wavelet Transform (DT-CWT) isused which provide perfect reconstruction over the traditional wavelet transform. Itemploys 2 real DWTs; the first DWT gives the real part of the transform while secondDWT gives the imaginary part. It is nearly shift invariant and directionally selective intwo and higher dimensions with limited redundancy. The DTCWT outperforms theDWT for applications like image denoising and enhancement. One of the advantagesof the DTCWT is that it can be used to implement 2D wavelet transforms that are moreselective with respect to orientation than is the 2D DWT. The 2D DTCWT producestwelve sub-bands at each scale, each of which are strongly oriented at distinct angles.Keywords: image ; denoising ; Dual Tree Complex wavelet; threshold前言近年来小波变换的快速发展，尤其是双树复数小波，很大程度上提升了图像处理算法的优化，双树复数小波最显著的优点是具有近似的平移不变性和更多的方向选择性。
双树小波变换与小波树稀疏联合的低场CS-MRI算法柴青焕;苏冠群;聂生东【摘要】压缩感知理论常用在磁共振快速成像上,仅采样少量的K空间数据即可重建出高质量的磁共振图像.压缩感知磁共振成像技术的原理是将磁共振图像重建问题建模成一个包含数据保真项、稀疏先验项和全变分项的线性组合最小化问题,显著减少磁共振扫描时间.稀疏表示是压缩感知理论的一个关键假设,重建结果很大程度上依赖于稀疏变换.本文将双树复小波变换和小波树稀疏联合作为压缩感知磁共振成像中的稀疏变换,提出了基于双树小波变换和小波树稀疏的压缩感知低场磁共振图像重建算法.实验表明,本文所提算法可以在某些磁共振图像客观评价指标中表现出一定的优势.【期刊名称】《波谱学杂志》【年(卷),期】2018(035)004【总页数】12页(P486-497)【关键词】低场磁共振成像;压缩感知;双树小波变换;小波树稀疏【作者】柴青焕;苏冠群;聂生东【作者单位】上海理工大学,医学影像工程研究所,上海 200093;上海理工大学,医学影像工程研究所,上海 200093;上海理工大学,医学影像工程研究所,上海 200093【正文语种】中文【中图分类】R445.2;O482.53磁共振成像（magnetic resonance imaging，MRI）技术是利用特定原子核（如氢核）的核磁共振（nuclear magnetic resonance，NMR）现象，通过定量检测物质中特定原子核共振过程中射频能量吸收和发射的规律，来显示物体内部结构及特征的一种成像技术．MRI扫描仪根据主磁场强度的不同，可分为低场、中场、高场、超高场等．一般情况下，低场≤0.3 T、0.3 T≤中场≤1.0 T、1.0 T≤高场≤3 T、超高场≥3 T．高场和超高场磁共振设备磁源多采用超导磁体产生高强磁场，采集到的信号幅值强、信噪比高．但超导磁体对设备环境要求高，且购买、安装和维护费用非常昂贵[1]．中场磁共振设备磁源采用超导磁体、永磁体或电磁，采集到的信号信噪比较高，一般用于医用MRI等．低场磁共振设备的磁源一般采用永磁体或电磁（室温下可使用），采集到的信号幅值很弱、信噪比极低．但是与中高场相比，低场磁共振设备具有便携、廉价、制造维护应用成本低，测量方便、快速、无损、无副作用等优点，已被广泛应用在食品化工、在线无损检测、聚合物材料、灾害防治、生物医药、能源勘探、石油测井、岩心分析等各种领域[2]，并逐渐成为一种必要的分析检查手段．近年来，随着电磁技术和计算机技术的发展，低场磁共振成像（low-field MRI，LF MRI）也迅速发展起来，高性能的低场开放型永磁MRI系统开始受到人们的青睐[3]．本文旨在研究低场磁共振图像重建算法，改善其图像质量、提高其成像速率．压缩感知（compressive sensing，CS）是一个新的信号采样和压缩理论．它利用信号在特定变换域上的稀疏性及可压缩性，能在远小于Nyquist采样率的条件下，从少量的信号中通过非线性重建算法恢复出原始信号[4]．传统的MRI技术需要很长的扫描时间，重建图像的时间较长．而将CS理论应用于MRI（CS-MRI）能够极大地减少傅里叶变换域的采样数据，降低扫描时间；且CS-MRI在不显著降低图像质量的情况下克服了MRI成像速度慢的缺点；验证了从少量稀疏测量数据重建出图像的可能性[5]．本文提出的低场CS-MRI重建技术借鉴了现已发展成熟的中高场CS-MRI技术，吸收其优点，并规避其缺点．对于CS-MRI，有两个关键点需要进一步研究．第一：稀疏变换，基于CS理论的MRI技术能从欠采样的稀疏图像数据中重建出原始图像，而磁共振图像本身不具有稀疏性，因此需要将图像在稀疏域进行稀疏表示．传统CS-MRI常用的稀疏基有小波变换、全变分（total variation，TV）等，然而它们有时不能提供足够的稀疏表示．如：小波变换（wavelet transform，WT）有以下缺点：具有转移敏感性、缺乏方向性、缺乏相位信息，且不能稀疏表示二维信号的一维奇异点，很难捕捉到规则的图像轮廓[6]．然而，通过改善可以弥补WT的缺点，如文献[7]中使用双树小波变换代替WT作为稀疏基，得到了好的稀疏效果．第二：重建算法，由于CS-MRI的重建过程是一个病态问题，一般使用正则约束项去求得病态问题的稳定解．例如：文献[8]提出利用共轭梯度算法去解TV和小波稀疏联合的正则化目标函数问题．文献[9]提出了一种用来解决与文献[5]同样问题的部分傅里叶系数重建（reconstruction partial Fourier，RecPF）算法．文献[10]提出一种将TV和CS应用于MRI的算法（TVCMRI），它采用一种算子分裂算法求解重建问题．文献[11]提出的快速复合分裂算法（fast composite splitting algorithm，FCSA）将正则化问题分裂为两个简单的子问题，每个子问题再通过快速收缩算法求解[11,12]，使正则化问题更有效的得到解决．最近，一种先进的小波树稀疏MRI重建算法（wavelet tree sparsity MRI，WaTMRI）[13]也已经被提出，它在TV和小波稀疏联合正则项基础上增加了小波树稀疏正则项．与先前的算法相比，WaTMRI算法使测量的复杂性由O(K+Klogn)减少到O(K+logn) [14]．其中K为信号稀疏度，n为原始信号的大小.基于以上理论，本文提出了一种双树小波变换与小波树稀疏联合的低场CS-MRI算法，称为DT-WaTMRI，在提高成像速度的同时重建出高分辨率的磁共振图像．以下部分是对双树小波变换与小波树稀疏联合的低场CS-MRI算法的详细阐述．CS的基本问题是从欠采样的线性观测向量中恢复出信号x，其中、、，（F是测量矩阵，R表示实数，n、m分别表示信号x、y的维度大小）；从数学角度求解以上欠定线性系统有无穷多解．然而，根据CS理论，在稀疏的假设下信号能通过(1)式所示的优化问题被恢复．其中是l0范数，表示向量中值不为0的个数．由于最小l0范数的优化问题是NP-hard问题，难以解决，Chrétien等人[15]证明可以利用凸的l1来近似非凸l0范数，如(2)式所示．(2)式是一个l1范数求解最小化的线性规划（linear programming, LP）问题，其中是变量x的绝对值之和，对其进行求解即可恢复出原始信号．CS-MRI就是将欠采样得到的少量K空间数据，进行某种稀疏变换，然后再结合(2)式的l1范数凸优化问题，最终重构出图像．CS-MRI的数学模型[8]可表示为(3)式其中x是待重建的磁共振图像；Y表示线性稀疏变换；表示欠采样傅里叶变换；y 是从MRI扫描仪欠采样测量得到的K空间数据；e控制重建数据和测量数据的保真度，代表测量数据的噪声水平[16]．利用拉格朗日乘子法将(3)式转换为无约束优化问题形式，如(4)式所示其中l为正则化参数，对(4)式进行求解即可重建出磁共振图像．近年来提出的各种CS-MRI方法都是基于改变(4)式的稀疏变换基或加入不同的正则项演化而来．目前，求解(4)式的CS-MRI算法多采用TV和小波稀疏约束项线性结合的方式来重建磁共振图像，如(5)式所示．其中a和b是两个正则参数；Y是小波变换；，和分别为图像水平方向和垂直方向的梯度算子，i、j表示磁共振图像在x、y方向像素位置．1.2.1 双树小波变换理论稀疏变换在CS理论中扮演着至关重要的角色，CS-MRI重建的结果很大程度上取决于稀疏变换．一般来说，信号越稀疏重建所需的测量值越少．除此之外，当存在测量不足和/或噪声时，重构图像中伪影的消除也依赖于稀疏转换[17]．双树小波变换（dual-tree complex wavelet transform, DT-CWT）是离散WT的增强形式，它能通过使用复值小波和尺度函数来弥补离散WT的缺点[6]. DT-CWT在二维或高维信号中具有移位不变性和定向选择性，能够很好地解决传统离散WT带来的在CS-MRI重建图像时产生边缘伪影、混淆现象及不具方向选择性等缺点．这使得DT-CWT成为CS-MRI中有效的稀疏变换工具．1.2.2 小波树稀疏理论自然数据（信号或图像）的小波系数的大部分系数值近似为0，只有少量的系数有较大值，具有很好的稀疏性．如：二维图像的小波系数服从四叉树结构，其在粗尺度上的系数被称为根节点，在细尺度上的系数被称为叶节点，每一根节点下面跟着四个细尺度上的子节点[18]，具有很高的稀疏性．图1所示为图像的小波四叉树结构．小波树结构为CS-MRI提供了良好的技术优势[19]．根据凸优化理论，结构稀疏可以通过凸优化算法求解；原始信号的先验信息，如组或图结构能被利用，还可以进一步降低采样率．当前提出的一些算法已经开始利用小波系数树结构来改进CS重构，这些算法很好的利用了树结构的理论，如：树结构的母-子关系，母系数有一个大的或小的值，那么它的子系数值往往与其保持一致．以上先验知识若能被完全利用，仅仅的观测数据就能恢复出满足树稀疏结构的原始数据，相比CS理论中标准K稀疏数据需要的采样数据才足够恢复出长度为n的原始数据，采样量减少了一半；当n无限大时，利用树结构的优势就更加显著了．文献[14]也验证了小波树结构能帮助提高成像速度．因此，本文提出了一种双树小波变换、小波树和梯度稀疏联合的新的稀疏变换模型，并将其应用于LF MRI，以便在提高其成像速度的同时，改善其图像质量．该模型中的稀疏变换双树小波变换和小波树稀疏被建模为一个重叠组正则项、TV作为惩罚项用来抑制噪声．TV惩罚项和l1范数约束项实际上是作为树形结构模型的补充，使新提出的稀疏变换模型对实际低场磁共振图像更具鲁棒性．实现以上条件的目标函数如(6)式所示：式中F表示双树复小波变换系数，G表示所有母-子系数组的集合，g为G的一个组，是属于母-子系数组的所有数据．现有的CS-MRI模型主要是求解最小二乘保真项、TV与小波稀疏正则项线性组合的目标函数最小化的问题．最早使用的是CG下降算法求解该问题．而TVCMRI和RecPF算法分别用作一种算子分裂算法和变量分裂算法来求解该问题；FCSA算法将原始问题分解为两个简单的子问题并用快速迭代收缩阈值算法（fast iterative shrinkage threshold algorithm, FISTA）分别求解它们．上述算法都是CS-MRI最常用的算法，但是它们都没有利用树状结构这一先验知识来增强其算法性能．文献[12]提出的WaTMRI算法能求解树结构、小波稀疏和TV联合的凸优化CS-MRI模型问题．并且该算法收敛速度快，每次迭代仅花费时间；而且在实际图像应用中，优于其它CS-MRI算法和几种一般的基于树的算法．因此本文也采用WaTMRI算法来求解双树小波变换、小波树和梯度稀疏联合的CS-MRI模型，称为DT-WaTMRI算法．DT-WaTMRI算法与先前CS-MRI的重建算法的主要不同是将稀疏变换基由小波变换改变为双树小波变换．这个算法不能用凸优化理论有效地求解．因此引入变量z约束x重叠结构，使(6)式变为非重叠凸优化问题，令，则(6)式衍变为下式：对于x子问题，我们可以联合(7)式中第一个和最后一个二次项惩罚项来求解，而z子问题则通过分组软阈值来解决，剩下部分与FCSA算法类似，可通过迭代方案有效解决．z子问题如(8)式所示：其中gi表示第i组，s表示全部组数．z子问题有一个封闭式的软阈值解决方案如(9)式所示：式中．令，这是一个有着莱布尼兹常数（Lipschitz, Lf）的凸平滑函数；再令，，和是凸的非平滑函数．为方便起见，将(9)式表示为(10)式：本文提出的DT-WaTMRI算法的主要步骤为：Input:For k=1 to N doEnd for算法中，N表示最大迭代次数，x表示待重建的图像，r迭代时存放相应值的中间变量，t表示适当的迭代步长；，表示函数f在点处的梯度，其中表示部分逆傅里叶转换．对任意标量，近端映射（proximal map，prox）归一化连续凸函数的定义如(11)式所示：(12)式将图像的像素值规范化到范围内，且，若不执行该操作，则重建图像由于存在负数像素数值而出现伪影．为验证所提出模型的良好性能，我们将DT-WaTMRI与FCSA和WaTMRI算法的重建结果进行比较．所有的算法都在处理器为Intel(R)Core(TM)i5-4460 CPU @3.20GHz的Dell Inspiron3847的Matlab2016b环境中运行．实验测试数据分为三组，分别为Matlab自带体模仿真数据，以及使用MesoMR23-040H-1的LF MRI扫描仪（苏州纽迈分析仪器有限公司）采集的哈密瓜及小鼠腹部的大小为256*128的全K空间数据．MRI扫描仪的磁场强度为0.5 T，脉冲序列选用自旋回波（spin echo，SE）序列，数据采集时实验参数分别设置如下：序列重复时间（repetition time, TR）为2 000 ms、回波时间（echo time, TE）为100 ms、扫描次数（number of scan, NS）为4．根据设定参数在LF MRI扫描仪上进行回波扫描得到所需数据．基于CS理论的MRI技术包括图像的稀疏表示、观测矩阵设计和算法重构三个方面．CS算法仅需采集Nyquist采样定理要求数据量的30％，甚至更少的数据重建原始信号．而磁共振图像的数据空间是K空间，K空间中存储在不同位置的磁共振信号对图像的对比度和空间分辨率的贡献有所差异：K空间的中心部分，信号强度大、空间频率低，形成图像对比度的权重大；K空间的外围部分，信号强度小、空间频率高，形成图像空间分辨率的权重大．因此CS-MRI对K空间数据欠采样选用变密度的随机测量矩阵，在K空间中心部分的采样权重多于K空间外围部分，即能加快采样速率又能保证成像质量．随机测量矩阵能很大概率的满足约束等距性质（restricted isometry property, RIP），RIP性质是CS理论完整恢复出原始信号的保证．但是随机测量矩阵硬件实现困难，所以当前实际应用中测量矩阵一般选用如图2所示（白色区域表示采样，黑色区域表示未采样）的伪随机确定性测量矩阵．由K空间的冗余特性采用图2所示的任意一种变密度的采样模式，都能缩短采样时间、提高成像速度、减少运动伪影．但是，不同采样脉冲序列产生不同的K空间模型，而不同的K空间模型需选用不同的随机变密度采样模式．由于实验所用采样脉冲序列为多回波SE序列，得到的K空间为笛卡尔单向平行轨迹，因此采用图2(a)所示的一维变密度采样方式．我们定义随机变密度欠采样矩阵的欠采样率为m/n（m和n分别为K空间的行、列大小），为方便起见，实验中所有算法的欠采样率均设为30%．利用DT-WaTMRI算法求解重建目标函数时，使其中的稀疏约束采用两级分解的Daubechies小波，正则化参数l、b、g分别设为0.1、0.1、0.025，最大迭代次数设置为200．在评价重建图像质量时，信噪比是一项重要的技术指标，但信噪比高的图像还不能确保将两个图像相邻区域或结构有效地区分开来，图像还必须在特定的组织和周围组织间有足够的对比度. 因此，本文对实验中不同算法的重建性能评估采用均方根误差（root mean square error, RMSE）、峰值信噪比（peak signal to noise ratio, PSNR）、方差（variance, Var）等指标来衡量.其公式分别如(13)、(14)、(15)式所示.公式(13)、(14)、(15)中表示全采样重建图像，表示欠采样重建图像，表示欠采样重建图像平均灰度值，M、N分别代表磁共振图像在x、y方向像素点个数，i、j 分别代表磁共振图像在x、y方向像素点位置.图3~5所示为不同数据的重建结果，图3(a)、4(a)、5(a)分别为不同样品的全采样重建图像、变密度采样矩阵、欠采样零填充重建图像；其中全采样图像作为欠采样重建图像的评价标准，变密度采样矩阵表示变密度欠采样K空间数据的采样方式，欠采样零填充重建图像是对欠采样的K空间数据不使用任何算法直接进行傅里叶变换得到的重建图像，零填充重建图像也作为不同重建算法欠采样重建图像的对比.图3(b)、4(b)、5(b)分别为使用FCSA、WaTMRI、DT-WaTMRI算法对体模数据、哈密瓜横断面数据、小鼠腹部横断面图像数据的重建结果. 图5中小鼠腹部横断面数据重建图像太小主要是因为对整个小鼠样本进行扫描时，受限于LF MRI扫描仪放置样品的口径，在前期处理中设定了成像视野（field of view, FOV），后续不易改动，但是不影响图像分辨率和后续的说明与结论.图6所示为使用FCSA、WaTMRI、DT-WaTMRI三种算法对体模数据、哈密瓜数据、小鼠腹部横断面数据的的重建结果的客观评价；图6(a)~(c)是对体模数据的重建结果的客观评价，图6(d)~(f)是对哈密瓜横断面数据的重建结果的客观评价，图6(g)~(i)是对小鼠腹部横断面图像数据的重建结果的客观评价；图6(a)、(d)、(g)为使用FCSA、WaTMRI、DT-WaTMRI三种算法欠采样重建图像与全采样重建图像的RMSE随迭代次数的变化情况，(b)、(e)、(h)为使用FCSA、WaTMRI、DT-WaTMRI三种算法欠采样重建图像与全采样重建图像的PSNR随迭代次数的变化情况，(c)、(f)、(i)为使用FCSA、WaTMRI、DT-WaTMRI三种算法欠采样重建图像的Var随迭代次数的变化情况.对CS-MRI来说重建时间也是评价算法性能的一项重要指标，因此本文也对不同算法的重建时间做了评价．表1列举了不同算法对不同样本的重建时间（单位为s）．图3、4、5所示为不同样本磁共振图像的重建结果，图3(a)、4(a)、5(a)中的零填充重建图像是不同样本的全采样K空间数据经一维变密度随机欠采样后得到欠采样的K空间数据，对欠采样的K空间中缺失的数据简单填零后直接傅里叶变换的重建结果，由结果可知零填充重建所得到的重建图像伪影比较严重．图3(b)、4(b)、5(b)是利用不同的CS-MRI重建算法对不同样本的欠采样数据进行重建的结果；由重建结果可知本文提出的DT-WaTMRI算法与FCSA算法和WaTMRI算法相比并没有显著改善图像的重建质量，而且使用本文算法对不同样本数据的重建结果还出现了类似伪影的竖条纹，尤其在体模数据中表现明显，这可能是因为体模数据是不含任何噪声的仿真数据，当用随机变密度矩阵对其进行欠采样，相当于原始无噪声数据矩阵与变密度矩阵相乘，造成了类似于条状伪影的数据缺失；也可能是本文算法还有需要改善的地方．但是本文算法与FCSA算法和WaTMRI算法相比，确实提高了重建图像的对比度和亮度，由图3也可以看出本文算法重建出图像的清晰度程度实际上是与原图最接近的，这点可由重建图像的评价指标Var验证.对重建出的结果图像可以通过视觉观察来判断其优劣，然而视觉观察往往带有选择性和主观性.为了提高说服力，图6采用客观评价指标如RMSE、PSNR、Var等，来判断不同算法的重建结果的优劣．RMSE是用来计算全采样重建图像和欠采样重建图像的均方值，然后通过均方值的大小来确定重建图像的失真程度.对不含噪声的体模数据来说，DT-WaTMRI算法的RMSE小于FCSA和WaTMRI算法；对含有噪声的哈密瓜和小鼠腹部低场磁共振图像来说，DT-WaTMRI算法的RMSE 略微小于FCSA和WaTMRI．PSNR是评价画质的客观测量法，值越大图像质量越好．对不含噪声的体模数据来说，DT-WaTMRI算法的PSNR大于FCSA和WaTMRI算法；对含有噪声的哈密瓜和小鼠腹部低场磁共振图像来说，DT-WaTMRI算法的PSNR略大于FCSA和WaTMRI算法．Var反应重建图像高频部分的大小，图像的对比度越大，方差就越大．对不含噪声的体模数据来说，DT-WaTMRI算法重建图像的Var大于FCSA和WaTMRI算法；对含有噪声的哈密瓜低场磁共振图像来说，DT-WaTMRI算法重建图像的Var值与FCSA算法相近；对含有噪声的小鼠腹部低场磁共振图像来说，DT-WaTMRI算法重建图像的Var 稍大于WaTMRI算法.表1中对所有实验样本来说，稀疏变换采用传统小波变换的FCSA算法重建出磁共振图像所用时间比稀疏变换采用传统小波变换与结构稀疏联合的WaTMRI算法少，这是由于算法中引入了小波树结构增加了计算冗余性．而稀疏变换采用双树小波与小波树结构联合的DT-WaTMRI算法重建出磁共振图像所用时间比WaTMRI 算法少，在保证重建质量的同时提高了成像效率．因此结合理论和实验结果可知本文提出的DT-WaTMRI算法改进了现有的重建算法，在一定程度上提高了磁共振图像的重建质量．在本文中，我们提出了一种基于双树小波变换与小波树结构联合的DT-WaTMRI 重建算法，并用WaTMRI算法加上一定的数据一致性约束对其进行求解.实验结果表明，该算法与其它算法在PSNR、RMSE和Var三个评价指标上，具有较好的性能；且与WaTMRI树稀疏算法相比，该算法重建计算消耗时间更少，适用于低场磁共振图像重建.【相关文献】[1] 赵刚. 低场磁共振分析仪的磁体和探头的设计[D]. 青岛: 山东科技大学, 2005.[2] WANG H M, NIE S D, WANG Y J. The research progress of de-noising methods in low-field NMR signal[J]. Chinese Journal ofMedica Physics, 2013, 30(4): 4261-4265. 王红敏, 聂生东, 王远军. 低场核磁共振信号降噪方法研究进展[J]. 中国医学物理学杂志, 2013, 30(4): 4261-4265.[3] 王鹤. 低场磁共振系统中若干技术问题的研究[D]. 上海: 华东师范大学, 2007.[4] SONG Y, XIE H B, YANG G. Segmentation dictionary learning algorithm for compressed sensing MRI[J]. Chinese J Magn Reson, 2016, 33(4): 559-569. 宋阳, 谢海滨, 杨光. 用于压缩感知磁共振成像的分割字典学习算法[J]. 波谱学杂志, 2016, 33(4): 559-569.[5] HUANG L J, SONG Y, ZHAO X C, et al. A new nombination scheme of GRAPPA and compressed sensing for accelerated magneticresonance imaging[J]. Chinese J Magn Reson, 2018, 35(1): 31-39.黄丽洁, 宋阳, 赵献策, 等. 一种结合并行成像和压缩感知的快速磁共振成像新方法[J]. 波谱学杂志, 2018, 35(1): 31-39.[6] SELESNICK I W, BARANIUK R G, KINGSBURY N C. The dual-tree complex wavelet transform[J]. IEEE Signal Proce Mag,2005, 22(6): 123-151.[7] KIM Y, ALTBACH M, TROUARD T, et al. Compressed sensing using dual-tree complex wavelet transform[C]. Proceedings of the International Society for Magnetic Resonance in Medicine, 2009.[8] LUSTIG M, DONOHO D, PAULY J M. Sparse MRI: The application of compressed sensing for rapid MR imaging[J]. Magn Reson Med, 2007, 58(6): 1182-1195.[9] YANG J F, ZHANG Y, YIN W T. A fast alternating direction method for TVL1-L2 signal reconstruction from partial fourier data[J]. IEEE J STSP, 2010, 4(2): 288-297.[10] MA S Q, YIN W T, ZHANG Y, et al. An efficient algorithm for compressed MR imaging using total variation and wavelets[C]. 2008 IEEE Conference on Computer Vision and Pattern Recognition, Anchorage: 2008.[11] HUANG J Z, ZHANG S T, METAXAS D. Efficient MR image reconstruction for compressed MR imaging[J]. Med Image Anal, 2011, 15(5): 670-679.[12] BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM J Imaging Sci, 2009, 2(1): 183-202.[13] CHEN C, HUANG J Z. Compressive sensing MRI with wavelet tree sparsity[C].International Conference on Neural Information Processing Systems, Lake Tahoe: 2012.[14] CHEN C, HUANG J Z. The benefit of tree sparsity in accelerated MRI[J]. Med Image Anal, 2014, 18(6): 834-842.[15] CHRÉTIEN S. An alternating l1 approach to the compressed sensing problem[J]. IEEE Signal Proc Let, 2010, 17(2): 181-184.[16] LUSTIG M, DONOHO D L, SANTOS J M, et al. Compressed sensing MRI[J]. IEEE Signal Proc Mag, 2008, 25(2): 72-82.[17] ZHU Z, WAHID K, BABYN P, et al. Compressed sensing-based MRI reconstruction using complex double-density dual-tree DWT[J]. Int J Biomed Imaging, 2013: 907501. [18] HUANG J Z, ZHANG T, METAXAS D. Learning with structured sparsity[J]. J Mach Learn Res, 2011, 12(7): 3371-3412.[19] BACH F, JENATTON R, MAIRAL J, et al. Structured sparsity through convex optimization[J]. Statist Sci, 2012, 27(4): 450-468.。
一、题目：双树复小波变换二、目的：双树复小波和实小波变换的比较三、算法及其实现：提取阶梯型边界点1．算法。
幅值:相位：2．代码实现。
我采用Matlab函数编程实现。
具体程序见shift_test_2D.m，drawcirc.m，setfig.m，dtwavexfm2.m，dtwaveifm2.m,waveifm2.m,wavexfm2.mSkelMap.m。
设和分别是双正交对偶尺度函数与对偶小波，，，和实部：虚部：双树复小波变换可以通过离散小波变换DWT实现：一个DWT产生实部，另一个产生虚部。
四、实现工具：Matlab五、程序代码：（1）shift_test_2D.m：% shift_test_2D.m%% M-file to perform a 4-level wavelet transform on a circle using Q-shift% dual wavelet tree and DWT, and to compare shift invariance properties.%% Nick Kingsbury, Cambridge University, May 2002.clear allclose all% Draw a circular disc.x = round((drawcirc(64,1,0,0,256) - 0.5) * 200);setfig(1);colormap(gray(256))image(min(max(x+128,1),256));set(gca,'position',[0.1 0.25 .25 .5]);()ctψ=,,j,c h gf fψψψ=+|(,)|cd j n=(,)arctan(,)gchd j nd j n∠= ⎪⎝⎭,h hφφ%,h hψψ%()h n()h n%()1h n()1h n%11()()(2)()()(2)()()(2)()()(2)h hnh hnh hnh hnt h n t nt h n t nt h n t nt h n t nφφψφφφψφ=-=-=-=-%%%%%%11()()(2)()()(2)()()(2)()()(2)g gng gng gng gnt g n t nt g n t nt g n t nt g n t nφφψφφφψφ=-=-=-=-%%%%%%axis('off');axis('image');% draw(xx);title('Input (256 x 256)','FontSize',14); drawnow% Do 4 levels of CWT.[Yl,Yh] = dtwavexfm2(x,4,'near_sym_b','qshift_b');% Loop to reconstruct output from coefs at each level in turn.% Starts with the finest level.titl = ['1st';'2nd';'3rd';'4th';'Low'];yy = zeros(size(x) .* [2 3]);yt1 = 1:size(x,1); yt2 = 1:size(x,2);for mlev = 1:5,mask = zeros(6,5);mask(:,mlev) = 1;z = dtwaveifm2(Yl*mask(1,5),Yh,'near_sym_b','qshift_b',mask);figure;draw(z);drawnowyy(yt1,yt2) = z;yt2 = yt2 + size(x,2)/2;end% disp('Press a key ...')% pause% Now do same with DWT.% Do 4 levels of Real DWT using 'antonini' (9,7)-tap filters. [Yl,Yh] = wavexfm2(x,4,'antonini');yt1 = [1:size(x,1)] + size(x,1); yt2 = 1:size(x,2);for mlev = 1:5,mask = zeros(3,5);mask(:,mlev) = 1;z = waveifm2(Yl*mask(1,5),Yh,'antonini',mask);figure;draw(z);drawnowyy(yt1,yt2) = z;yt2 = yt2 + size(x,2)/2;endfigure;setfig(gcf);colormap(gray(256))image(min(max(yy+128,1),256));set(gca,'position',[0.1 0.1 .8 .8]);axis('off');axis('image');hold onplot(128*[[1;1]*[1:4] [0;6]]+1,128*[[0;4]*[1 1 1 1] [2;2]]+1,'-k');hold offtitle('Components of reconstructed ''disc'' images','FontSize',14);text(-0.01*size(yy,2),0.25*size(yy,1),'DT CWT','horiz','r');text(0.02*size(yy,2),1.02*size(yy,1),'wavelets:','horiz','r','vert','t');text(-0.01*size(yy,2),0.75*size(yy,1),'DWT','horiz','r');for k=1:4, text(k*128-63,size(yy,1)*1.02,sprintf('level %d',k),'FontSize',14,'horiz','c','vert','t'); endtext(5*128+1,size(yy,1)*1.02,'level 4 scaling fn.','FontSize',14,'horiz','c','vert','t');drawnow% print -deps circrecq.epsdisp('Press a key to see perfect reconstruction property ...')pause% Accumulate the images from lowband upwards to show perfect reconstruction.sy = size(x,2)/2;for mlev = 4:-1:1,yt2 = [1:sy] + (mlev-1)*sy;yy(:,yt2) = yy(:,yt2) + yy(:,yt2+sy);endfigure;setfig(gcf);colormap(gray(256))image(min(max(yy+128,1),256));set(gca,'position',[0.1 0.1 .8 .8]);axis('off');axis('image');title('Accumulated reconstructions from each level of DT CWT ','FontSize',14);text(size(yy,2)*0.5,size(yy,1)*1.02,'Accumulated reconstructions from each level of DWT ','FontSize',14,'hor','c','vert','t');drawnowreturn（2）function p = drawcirc(r,w,dx,dy,N)% function p = drawcirc(r,w,dx,dy,N)% Generate an image of size N*N pels, containing a circle% radius r pels and centred at dx,dy relative% to the centre of the image. The edge of the circle is a cosine shaped% edge of width w (from 10 to 90% points).x = ones(N,1) * (([1:N] - (N+1)/2 - dx)/r);y = (([1:N]' - (N+1)/2 - dy)/r) * ones(1,N);p = 0.5 + 0.5 * sin(min(max((exp(-0.5 * (x + y )) - exp(-0.5))*(r*3/w), -pi/2), pi/2));return（3）function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree complex wavelet (DTCWT)% 2-D reconstruction.%% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);%% Yl -> The real lowpass image from the final level% Yh -> A cell array containing the 6 complex highpass subimages for each level.%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,% (only 6,6 non-zero taps).% 'qshift_a' => Q-shift 10,10 tap filters,% (with 10,10 non-zero taps, unlike qshift_06).% 'qshift_b' => Q-Shift 14,14 tap filters.% 'qshift_c' => Q-Shift 16,16 tap filters.% 'qshift_d' => Q-Shift 18,18 tap filters.%% gain_mask -> Gain to be applied to each subband.% gain_mask(d,l) is gain for subband with direction d at level l.% If gain_mask(d,l) == 0, no computation is performed for band (d,l).% Default gain_mask = ones(6,length(Yh)).%% Z -> Reconstructed real image matrix%%% For example: Z = dtwaveifm2(Yl,Yh,'near_sym_b','qshift_b');% performs a 3-level reconstruction from Yl,Yh using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, May 2002a = length(Yh); % No of levels.if nargin < 5, gain_mask = ones(6,a); end % Default gain_mask.if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_exist = exist([biort '.mat']);qshift_exist = exist([qshift '.mat']);if biort_exist == 2 & qshift_exist == 2; %Check to see if the inputs exist as .mat files load (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEIFM2.');endcurrent_level = a;Z = Yl;while current_level >= 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do even Qshift filters on columns.y1 = colifilt(Z,g0b,g0a) + colifilt(lh,g1b,g1a);y2 = colifilt(hl,g0b,g0a) + colifilt(hh,g1b,g1a);% Do even Qshift filters on rows.Z = (colifilt(y1.',g0b,g0a) + colifilt(y2.',g1b,g1a)).';% Check size of Z and crop as required[row_size col_size] = size(Z);S = 2*size(Yh{current_level-1});if row_size ~= S(1) %check to see if this result needs to be cropped for the rows Z = Z(2:row_size-1,:);endif col_size ~= S(2) %check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1);endif any(size(Z) ~= S(1:2)),error('Sizes of subbands are not valid for DTW A VEIFM2');endcurrent_level = current_level - 1;endif current_level == 1;lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do odd top-level filters on columns.y1 = colfilter(Z,g0o) + colfilter(lh,g1o);y2 = colfilter(hl,g0o) + colfilter(hh,g1o);% Do odd top-level filters on rows.Z = (colfilter(y1.',g0o) + colfilter(y2.',g1o)).';endreturn%==================================================================== ======================% ********** INTERNAL FUNCTION **********%==================================================================== ======================function x = c2q(w,gain)% function z = c2q(w,gain)% Scale by gain and convert from complex w(:,:,1:2) to real quad-numbers in z. %% Arrange pixels from the real and imag parts of the 2 subbands% into 4 separate subimages .% A----B Re Im of w(:,:,1)% | |% | |% C----D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2));if any(w(:)) & any(gain)sc = sqrt(0.5) * gain;P = w(:,:,1)*sc(1) + w(:,:,2)*sc(2);Q = w(:,:,1)*sc(1) - w(:,:,2)*sc(2);t1 = 1:2:size(x,1);t2 = 1:2:size(x,2);% Recover each of the 4 corners of the quads.x(t1,t2) = real(P); % a = (A+C)*sc;x(t1,t2+1) = imag(P); % b = (B+D)*sc;x(t1+1,t2) = imag(Q); % c = (B-D)*sc;x(t1+1,t2+1) = -real(Q); % d = (C-A)*sc;endreturn（3）function [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort,qshift);% Function to perform a n-level DTCWT-2D decompostion on a 2D matrix X %% [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort,qshift);%% X -> 2D real matrix/Image%% nlevels -> No. of levels of wavelet decomposition%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,% (only 6,6 non-zero taps).% 'qshift_a' => Q-shift 10,10 tap filters,% (with 10,10 non-zero taps, unlike qshift_06).% 'qshift_b' => Q-Shift 14,14 tap filters.% 'qshift_c' => Q-Shift 16,16 tap filters.% 'qshift_d' => Q-Shift 18,18 tap filters.%%% Yl -> The real lowpass image from the final level% Yh -> A cell array containing the 6 complex highpass subimages for each level.% Yscale -> This is an OPTIONAL output argument, that is a cell array containing% real lowpass coefficients for every scale.%%% Example: [Yl,Yh] = dtwavexfm2(X,3,'near_sym_b','qshift_b');% performs a 3-level transform on the real image X using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, Sept 2001if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_exist = exist([biort '.mat']);qshift_exist = exist([qshift '.mat']);if biort_exist == 2 & qshift_exist == 2; %Check to see if the inputs exist as .mat filesload (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEXFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEXFM2.');endorginal_size = size(X);if ndims(X) >= 3;error(sprintf('The entered image is %dx%dx%d, please enter each image sliceseparately.',orginal_size(1),orginal_size(2),orginal_size(3)));end% The next few lines of code check to see if the image is odd in size, if so an extra ...% row/column will be added to the bottom/right of the imageinitial_row_extend = 0; %initialiseinitial_col_extend = 0;if any(rem(orginal_size(1),2)), %if sx(1) is not divisable by 2 then we need to extend X by adding a row at the bottomX = [X; X(end,:)]; %Any further extension will be done in due course.initial_row_extend = 1;endif any(rem(orginal_size(2),2)), %if sx(2) is not divisable by 2 then we need to extend X by adding a col to the leftX = [X X(:,end)]; %Any further extension will be done in due course.initial_col_extend = 1;endextended_size = size(X);if nlevels == 0, return; end%initialiseYh=cell(nlevels,1);if nargout == 3Yscale=cell(nlevels,1); %this is only required if the user specifies a third output component.endS = [];sx = size(X);if nlevels >= 1,% Do odd top-level filters on cols.Lo = colfilter(X,h0o).';Hi = colfilter(X,h1o).';% Do odd top-level filters on rows.LoLo = colfilter(Lo,h0o).'; % LoLoYh{1} = zeros([size(LoLo)/2 6]);Yh{1}(:,:,[1 6]) = q2c(colfilter(Hi,h0o).'); % Horizontal pairYh{1}(:,:,[3 4]) = q2c(colfilter(Lo,h1o).'); % Vertical pairYh{1}(:,:,[2 5]) = q2c(colfilter(Hi,h1o).'); % Diagonal pairS = [ size(LoLo) ;S];if nargout == 3Yscale{1} = LoLo;endendif nlevels >= 2;for level = 2:nlevels;[row_size col_size] = size(LoLo);if any(rem(row_size,4)), % Extend by 2 rows if no. of rows of LoLo are divisable by 4;LoLo = [LoLo(1,:); LoLo; LoLo(end,:)];endif any(rem(col_size,4)), % Extend by 2 cols if no. of cols of LoLo are divisable by 4;LoLo = [LoLo(:,1) LoLo LoLo(:,end)];end% Do even Qshift filters on rows.Lo = coldfilt(LoLo,h0b,h0a).';Hi = coldfilt(LoLo,h1b,h1a).';% Do even Qshift filters on columns.LoLo = coldfilt(Lo,h0b,h0a).'; %LoLoYh{level} = zeros([size(LoLo)/2 6]);Yh{level}(:,:,[1 6]) = q2c(coldfilt(Hi,h0b,h0a).'); % HorizontalYh{level}(:,:,[3 4]) = q2c(coldfilt(Lo,h1b,h1a).'); % VerticalYh{level}(:,:,[2 5]) = q2c(coldfilt(Hi,h1b,h1a).'); % DiagonalS = [ size(LoLo) ;S];if nargout == 3Yscale{level} = LoLo;endendendYl = LoLo;if initial_row_extend == 1 & initial_col_extend == 1;warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r The bottom row and rightmost column have been duplicated, prior to decomposition. \r\r ',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2)));endif initial_row_extend == 1 ;warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r Row number %d has been duplicated, and added to the bottom of the image, prior to decomposition. \r\r',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(1)));endif initial_col_extend == 1;warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r Col number %d has been duplicated, and added to the right of the image, prior to decomposition. \r\r',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(2)));endreturn%==================================================================== ======================% ********** INTERNAL FUNCTION **********%==================================================================== ======================function z = q2c(y)% function z = q2c(y)% Convert from quads in y to complex numbers in z.sy = size(y);t1 = 1:2:sy(1); t2 = 1:2:sy(2);j2 = sqrt([0.5 -0.5]);% Arrange pixels from the corners of the quads into% 2 subimages of alternate real and imag pixels.% a----b% | |% | |% c----d% Combine (a,b) and (d,c) to form two complex subimages.p = y(t1,t2)*j2(1) + y(t1,t2+1)*j2(2); % p = (a + jb) / sqrt(2)q = y(t1+1,t2+1)*j2(1) - y(t1+1,t2)*j2(2); % q = (d - jc) / sqrt(2)% Form the 2 subbands in z.z = cat(3,p-q,p+q);return（4）function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree complex wavelet (DTCWT)% 2-D reconstruction.%% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);%% Yl -> The real lowpass image from the final level% Yh -> A cell array containing the 6 complex highpass subimages for each level.%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters, % (only 6,6 non-zero taps).% 'qshift_a' => Q-shift 10,10 tap filters,% (with 10,10 non-zero taps, unlike qshift_06).% 'qshift_b' => Q-Shift 14,14 tap filters.% 'qshift_c' => Q-Shift 16,16 tap filters.% 'qshift_d' => Q-Shift 18,18 tap filters.%% gain_mask -> Gain to be applied to each subband.% gain_mask(d,l) is gain for subband with direction d at level l.% If gain_mask(d,l) == 0, no computation is performed for band (d,l). % Default gain_mask = ones(6,length(Yh)).%% Z -> Reconstructed real image matrix%%% For example: Z = dtwaveifm2(Yl,Yh,'near_sym_b','qshift_b');% performs a 3-level reconstruction from Yl,Yh using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, May 2002a = length(Yh); % No of levels.if nargin < 5, gain_mask = ones(6,a); end % Default gain_mask.if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_exist = exist([biort '.mat']);qshift_exist = exist([qshift '.mat']);if biort_exist == 2 & qshift_exist == 2; %Check to see if the inputs exist as .mat files load (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEIFM2.');endcurrent_level = a;Z = Yl;while current_level >= 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do even Qshift filters on columns.y1 = colifilt(Z,g0b,g0a) + colifilt(lh,g1b,g1a);y2 = colifilt(hl,g0b,g0a) + colifilt(hh,g1b,g1a);% Do even Qshift filters on rows.Z = (colifilt(y1.',g0b,g0a) + colifilt(y2.',g1b,g1a)).';% Check size of Z and crop as required[row_size col_size] = size(Z);S = 2*size(Yh{current_level-1});if row_size ~= S(1) %check to see if this result needs to be cropped for the rows Z = Z(2:row_size-1,:);endif col_size ~= S(2) %check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1);endif any(size(Z) ~= S(1:2)),error('Sizes of subbands are not valid for DTW A VEIFM2');endcurrent_level = current_level - 1;endif current_level == 1;lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do odd top-level filters on columns.y1 = colfilter(Z,g0o) + colfilter(lh,g1o);y2 = colfilter(hl,g0o) + colfilter(hh,g1o);% Do odd top-level filters on rows.Z = (colfilter(y1.',g0o) + colfilter(y2.',g1o)).';endreturn%==================================================================== ======================% ********** INTERNAL FUNCTION **********%==================================================================== ======================function x = c2q(w,gain)% function z = c2q(w,gain)% Scale by gain and convert from complex w(:,:,1:2) to real quad-numbers in z.%% Arrange pixels from the real and imag parts of the 2 subbands% into 4 separate subimages .% A----B Re Im of w(:,:,1)% | |% | |% C----D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2));if any(w(:)) & any(gain)sc = sqrt(0.5) * gain;P = w(:,:,1)*sc(1) + w(:,:,2)*sc(2);Q = w(:,:,1)*sc(1) - w(:,:,2)*sc(2);t1 = 1:2:size(x,1);t2 = 1:2:size(x,2);% Recover each of the 4 corners of the quads.x(t1,t2) = real(P); % a = (A+C)*sc;x(t1,t2+1) = imag(P); % b = (B+D)*sc;x(t1+1,t2) = imag(Q); % c = (B-D)*sc;x(t1+1,t2+1) = -real(Q); % d = (C-A)*sc;endreturn(5)function [Yl,Yh,Yscale] = wavexfm2(X,nlevels,biort);% Function to perform a n-level DWT-2D decompostion on a 2-D matrix X.%% [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort);%% X -> real 1-D signal column vector (or matrix of vectors)%% nlevels -> No. of levels of wavelet decomposition%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% Yl -> The lowpass subband from the final level.% Yh -> A cell array containing the highpass subband for each level.% Yscale -> This is an OPTIONAL output argument, that is a cell array containing% the lowpass coefficients at every scale.%%% Example: [Yl,Yh] = wavexfm2(X,4,'near_sym_b');% performs a 4-level 2-D DWT on the real image X using the 13,19-tap filters.%% Nick Kingsbury, Cambridge University, May 2002if isstr(biort) % Check if the biort input is a stringbiort_exist = exist([biort '.mat']);if biort_exist == 2, % Check to see if the filter exists as a .mat fileload (biort);elseerror('Please enter the correct name of the Biorthogonal Filter, see help W A VEXFM2 for details.');endelseerror('Please enter the name of the Biorthogonal Filter as shown in help W A VEXFM2.');endL = size(X);if any(rem(L,2)), % ensure that X is an even length, thus enabling it to be extended if needs be.error('Size of X must be a multiple of 2');end%initialiseYh=cell(nlevels,1);if nargout == 3Yscale=cell(nlevels,1); % This is only required if the user specifies a third output component.endLoLo = X;for level = 1:nlevels;if rem(size(LoLo,1),4), % Check to see if height of LoLo is divisable by 4, if not extend.LoLo = [LoLo(1,:); LoLo; LoLo(end,:)];endif rem(size(LoLo,2),4), % Check to see if height of LoLo is divisable by 4, if not extend.LoLo = [LoLo(:,1) LoLo LoLo(:,end)];end% Do filters on rows.Lo = coldwtfilt(LoLo,h0o,0).';Hi = coldwtfilt(LoLo,h1o,1).';% Do filters on columns.LoLo = coldwtfilt(Lo,h0o,0).'; %LoLoYh{level} = zeros([size(LoLo) 3]);Yh{level}(:,:,1) = coldwtfilt(Hi,h0o,0).'; % HorizontalYh{level}(:,:,3) = coldwtfilt(Lo,h1o,1).'; % VerticalYh{level}(:,:,2) = coldwtfilt(Hi,h1o,1).'; % Diagonalif nargout == 3Yscale{level} = LoLo;endendYl = LoLo;return(6)function Z = waveifm2(Yl,Yh,biort,gain_mask);% Function to perform an n-level DWT 2-D reconstruction.%% Z = waveifm2(Yl,Yh,biort,gain_mask);%% Yl -> The real lowpass subband from the final level% Yh -> A cell array containing the complex highpass subband for each level.%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% gain_mask -> Gain to be applied to each subband.% gain_mask(l) is gain for wavelet subband at level l.% If gain_mask(l) == 0, no computation is performed for band (l).% Default gain_mask = ones(1,length(Yh)).%% Z -> Reconstructed real signal vector (or matrix).%%% For example: Z = waveifm2(Yl,Yh,'near_sym_b');% performs a reconstruction from Yl,Yh using the 13,19-tap filters.%% Nick Kingsbury, Cambridge University, May 2002nlevels = length(Yh); % No of levels.if nargin < 4, gain_mask = ones(3,nlevels); end % Default gain_mask.if isstr(biort) % Check if the biort input is a stringbiort_exist = exist([biort '.mat']);if biort_exist == 2, % Check to see if the filter exists as a .mat fileload (biort);elseerror('Please enter the correct name of the Biorthogonal Filter, see help W A VEIFM2 for details.');endelseerror('Please enter the name of the Biorthogonal Filter as shown in help W A VEIFM2.');endLoLo = Yl;for level = nlevels:-1:1, % Reconstruct levels in reverse order.if size(LoLo,1) ~= size(Yh{level},1) % If LoLo is not the same height as the next Yh => t1 was extended.LoLo = LoLo(2:size(LoLo,1)-1,:); % Therefore we have to clip LoLo so it is the same height as the next Yh.endif size(LoLo,2) ~= size(Yh{level},2) % If LoLo is not the same width as the next Yh => t1was extended.LoLo = LoLo(:,2:size(LoLo,2)-1); % Therefore we have to clip LoLo so it is the same width as the next Yh.endif any([size(LoLo) 3] ~= size(Yh{level})),error('Yh sizes are not valid for WA VEIFM2');endlh = Yh{level}(:,:,1) * gain_mask(1,level);hl = Yh{level}(:,:,3) * gain_mask(3,level);hh = Yh{level}(:,:,2) * gain_mask(2,level);% Do even Qshift filters on columns.y1 = coliwtfilt(LoLo,2*g0o,0) + coliwtfilt(lh,2*g1o,1);y2 = coliwtfilt(hl,2*g0o,0) + coliwtfilt(hh,2*g1o,1);% Do even Qshift filters on rows.LoLo = (coliwtfilt(y1.',2*g0o,0) + coliwtfilt(y2.',2*g1o,1)).';endZ = LoLo;return(7)function setfig(f)% function setfig(f)% Set figure dimensions for correct display of text to agree with% print -deps.figure(f)set(gcf,'position',[220 40 790 526],'DefaultTextFontSize',14);return六、运行结果：七、结果分析：。
