基于快速非局部均值和TV-L1模型的剪切波变换医学CT图像去噪方法与流程

本发明涉及于医学图像去噪领域，特别是涉及医学ct图像，具体涉及一种适用于医学ct图像的基于快速非局部均值和tv-l1模型的剪切波变换医学ct图像去噪方法。

背景技术：

随着科技的发展，在医学成像领域，超声成像、ct、mri等成像技术已应用于医学临床诊断中。计算机断层成像(computedtomography，又称为“电脑断层扫描”，简称ct)，是一种影像诊断学的检查。这一技术曾被称为计算机轴向断层成像(computedaxialtomography)。计算机断层扫描，利用计算机处理许多的组合x射线测量了从不同角度产生的横断面的特定区域扫描对象，允许用户看到对象的内部并且没有削减。由于ct成像技术检查为横断面成像，可以通过图像重建，任意方位显示组织或器官，对病变显示更为全面，防止遗漏；具有高密度分辨率，对有密度改变的细微病变也能显示出来，可以明确病变的性质；此外，ct具有无创、成像快等优势已经成为一种广泛使用且高度安全的医疗诊断技术。

均匀物体的影像中各象素的ct值参差不齐，图像呈颗粒性，影响密度分辨力，这种现象称ct的噪声。其来源有探测器方面的，如探测器的灵敏度，象素大小，层厚及x线量等。还有电子线路及机械方面的，重建方法及散乱射线等也会引起噪声。噪声与图像的质量成反比，因此要了解噪声产生的机制，尽量加以抑制。

本发明使用医学ct为研究对象.由于ct成像不免受到各种物理因素的影响，斑点噪声的存在严重影响了ct图像的质量，导致了医学图像质量较差。斑点噪声在图像上表现为空间域内相关的形状各异的小斑点，它将掩盖那些灰度差别很小的图像特征。对于临床医生而言，斑点噪声对他们的准确诊断造成了很大的干扰，特别是对于经验不是很丰富的医生造成的影响更大。因此，从临床应用的角度出发，需要研究对ct医学图像去噪方法，为医生做出更准确的诊断提供技术支持，降低人工诊断的风险。

综上所述，研究医学ct图像去噪方法具有非常重要的意义。

技术实现要素：

本发明要克服现有技术的小波分析在处理高维数据稀疏能力方面不足的缺点，提供了一种离散离剪切波医学ct图像去噪算法,用于降低医学ct图像上的无效噪声。在传统的技术中，小波变换的优点是可以非常有效的应用在图像去噪上而且很好的处理一维奇点，缺点是无法应对曲线和直线的突变。在这一点脊波变换能够有效的应用于捕捉线的奇异性，克服小波变换在这一方面的缺点，不过依旧无法非常有效的捕捉曲线的奇异性。从近几年的研究现状上来说，离散剪切波算法在医学图像上的应用，使医学图像去噪技术领域产生了一定程度的进步。shearlet变换具有多分辨率、方向性、局部性、各向异性，是图像最稀疏的表示，并且已经在图像处理领域有了广泛的应用。本发明正是将离散剪切波变换应用到医学ct图像的去噪当中，另外本发明还将tv-l1模型应用到图像分解当中，用分解得到的由灰度和边界组成的平滑变化的成分和剪切波分解得到的低频部分混合得到新的低频部分，而后进行后续的处理。其中l1特别适合于平滑成分和纹理成分的分解。本发明还应用了快速非局部均值去噪的方法发明了具有速度快、去噪明显的剪切波医学ct图像去噪声方法，最后通过仿真验证了方法的可行性与优化的效果。

与现有技术相比，本发明的创造性与优点在于：本发明将tv-l1模型用于去噪图像的分解，将分解得出的平滑成分与剪切波分解得出的低频成分进行了结合得到了新的低频成分；提出了快速的非局部均值去噪方法应用于降低图像无效噪声；在此基础上结合剪切波的特性克服了小波分析处理高维数据稀疏能力的不足及脊波不能很好捕捉曲线的奇异性不足。将此具有多分辨、多尺度、多方向性的方法应用于ct图像去噪可以更好的保护图像边缘信息让医师作出更好的诊断。

为使本发明的目的、技术方案和优点更加清晰，下面就对本发明的技术方案作进一步描述，基于快速非局部均值和tv-l1模型的剪切波变换医学ct图像去噪方法，包括以下六个步骤。

步骤1)建立医学ct图像的模型；

计算机断层扫描，利用x射线从多个不同的角度和方位扫描人体部位再由计算机处理不同的横断面得到重建图像，使用户看到特定区域的扫描对象。由于低强度的发射电流会产生足以影响观察判断的高斯噪声，降低生成图像的图像质量。

ct图像的模型主要由两部分组成，即有效的人体组织反射信号和无效的噪声信号，而噪声信号则包括乘性噪声和加性噪声，其中加性噪声相对于乘性噪声来说对ct图像的影响非常微小，在通常的处理中我们忽略加性噪声。这样一来ct电信号的通用模型可以表示为：

o(x,y)＝p(x,y)q(x,y)(1)

其中，x表示图像的横坐标，y表示图像的纵坐标，r(x,y)表示无噪声信号，n(x,y)表示相乘性噪声。

为了便于噪声的分离，需要对以上的模型进行对数变换将相乘的式(1)模型变换成相加模型：

log(o(x,y))＝log(p(x,y))+log(q(x,y))(2)

步骤2)对图像进行剪切波变换多尺度和多方向分解；

对上一步对数变换后便于噪声分离的ct图像应用多尺度分解，经分解后得到一张和原图等大的低频ct图像分量和多张同样与原图等大的高频ct图像分量。其中低频部分不做处理，在步骤3中会进行下一步处理，对于高频部分采用剪切滤波器作用到各个ct子带，对各尺度子带图像用剪切滤波器组进行方向分解。

步骤3)使用tv-l1模型分解图像为cartoon和texture部分；

运用tv-l1模型可以将图像分解成平滑的cartoon部分和保留纹理和噪声的texture部分，tv-l1模型可表示为：

其中u和f分别是输出图像和输入图像，tv-l1模型在不同的尺度下得到的结构不同，在尺度空间上能够有效的把不同的结构区分开。参数λ取值的大小对图像分解平滑有很大的作用，λ取值越小则平滑的越细致，cartoon部分的平滑效果非常干净。并且分解的留一部分，也就是细节纹理部分(texture)也非常清晰。这里取上一步中的低频部分和这里的cartoon部分相混合组成新的低频部分，为下一步的处理做好准备。

步骤4)对混合图应用快速非局部均值去噪；

对混合图应用利用积分图像技术加速得到的快速非局部均值去噪，得到新的低频子带。非局部均值充分利用了图像中的冗余信息，在尽可能的保留图像的细节特征的同时完成图像的去噪。非局部均值的基本思想是：当前像素的估计值的根据图像中和它有类似邻域结构的像素加权平均得来的。

设含噪声图像为v，去噪后的图像为u。u中像素点x处的灰度大小的获取可以表示为：

其中权值w(x,y)表示像素点x和y之间的相似度，由以x，y为中心的矩形邻域v(x)、v(y)间的距离||v(x)-v(y)||²决定：

z(x)是归一化系数，h是平滑参数，h越大去噪程度越深，但如果过大则会导致图像变得模糊；如果h越小，越能够保留足够多的边缘细节成分，但h太小又会有过多噪声。

由于原始非局部均值算法复杂度太高，计算耗时非常大，有必要对其进行改进。

这里采用快速非局部均值方法主要由积分图像技术加速，积分图像技术是一种计算图像中矩形区域和的快速方法，积分图像技术的长处只要先计算出积分图像我们就能够在常量时间范围内计算图像任何大小的矩形的和。在积分图像(integralimage-ii)上任意一个位置(x,y)上的ii(x,y)表示该点左上角矩形范围内所有像素的和，表示如下：

其中i(x,y)表示输入图像上相关位置的像素点。若给定图像i，则从上到下、从左到右计算得到和的积分图像公式表示如下：

ii(x,y)＝ii(x-1,y)+i(x,y)+ii(x,y-1)-ii(x-1,y-1)(8)

其中(x＜0||y＜0)时ii(x,y)＝0,i(x,y)＝0。得到积分图像之后，图像中任意矩形区域和通过如下公式计算：

sum(m,n)＝ii(x,y)+ii(u,v)-ii(x,v)-ii(u,y)(9)

其中矩形大小为：m＝x-u,n＝y-v。

在利用积分图像技术为非局部均值去噪加速的应用中，假如待处理图像一共有像素点n个，设置搜索窗口为d×d(d＝2×ds+1)大小，设置邻域窗口为d×d(d＝2×ds+1)大小，那么两个矩形邻域之间的计算相似度时间就是o(d²)，其中每一个像素点都需要计算与搜索窗口内任何一个像素之间的相似度，而搜索窗口内具有d²个像素，所以可得非局部均值nl-means的计算复杂度就是o(nd²d²)。想要提高nl-means的运行计算速度就好从其计算复杂度o(nd²d²)入手，就该计算复杂度而言可以从o(d²)入手，也就是从降低邻域间相似度的计算入手。在原始计算过程中，想要计算邻域之间的距离的时候需要对每一个像素点求得差值，遍历所求距离的两个邻域。首先构造一个关于像素差值的积分图像可以解决这个问题：

其中st(x)＝||v(x)-v(x+t)||²，由此在计算两个邻域v(x)和v(y)(y＝x+t)之间的距离的时候，就可以在常量时间内完成计算：

为了提高运算速度，降低空间复杂度，以上所讲的算法把偏移量当做最外层的循环，也就是说每次只需要在一个偏移方向上求取积分图像，并对该积分图像进行处理即可。而不需要一次性求取出所有积分图像。由此，原计算复杂度o(nd²d²)降低为o(nd²)，加快了算法运行速度。

步骤5)对高频子带的剪切波系数进行阈值收缩处理；

在对高频部分进行阈值收缩处理的时候，其阈值函数的选取在很大程度上决定这去噪效果的好坏。若阈值选取过小，虽然有用信息不回被去除但同时也会有过多噪声被保留下来，若阈值确定的太大，虽然噪声是被去除了，但是许多有效信息也被去除了。在传统滤波阈值函数算法有软阈值和硬阈值算法，但是在多尺度多方向的剪切波系数处理中无法得到非常有效的应用效果。对此，本发明采用新型的自适应阈值收缩算法。对得到的多个高频ct子带的剪切波系数采取自适应阈值收缩处理，本发明使用的针对ct医学图像的阈值函数：

式中，σn是噪声的标准差，m是对应变换域内变换系数的总体个数，tj代表j层的自适应参数，经剪切波分解后不同分解层的系数会有不同的分布，所以tj基于不同的j层及具体实验目标选取；在本发明中，512*512的ct图分解为一个低频部分和4个高频部分，阈值t≈5σn，对于细节较多的分解子带，可使用t≈4σn或t≈3σn。

步骤6)对处理后的系数进行剪切波逆变换；

经过阈值收缩以及快速非局部均值处理有就得到了去噪后的剪切波系数，想要获得去噪之后的ct图像，就需要对处理后的系数进行剪切波逆变换，由此得到方便医师诊断分析的清晰的ct图像。

本发明具有以下优点：

1.本发明将tv-l1模型用于去噪图像的分解，将分解得出的平滑成分与剪切波分解得出的低频成分进行了结合得到了新的低频成分；

2.本发明采用积分图像技术加快了非局部均值去噪方法，并将只应用于低频混合图像降低图像无用噪声，去噪声的同时保留边缘有效信息；

3.本发明采用自适应的剪切波系数阈值收缩算法，可以有效处理分解后的子带噪声；

4.本发明步骤明确结构简单，拥有完善的理论支持。

附图说明

图1为lena图片经tv-l1模型分解所得的cartoon部分和经分解处理所得的低频部分的混合图；

图2为本发明整体步骤流程图；

图3为案例分析整体流程；

图4为各种算法应用在lena图(σ＝40)的实验效果图以及原图和噪声图，图4a为原图，图4b是噪声图，图4c为fdct算法效果图，图4d是ffst算法效果图，图4e是本发明算法效果图；

图5为各种算法应用在头部ct图(σ＝40)的实验效果图以及原图和噪声图，图5a是噪声图，图5b为fdct算法效果图，图5c是ffst算法效果图，图5d是本发明算法效果图。

具体实施方式

以下结合附图对本发明做进一步说明。

本发明的基于快速非局部均值和tv-l1模型的剪切波变换医学ct图像去噪方法，包括以下步骤：

步骤1)建立医学ct图像的模型；

计算机断层扫描，利用x射线从多个不同的角度和方位扫描人体部位再由计算机处理不同的横断面得到重建图像，使用户看到特定区域的扫描对象。由于低强度的发射电流会产生足以影响观察判断的高斯噪声，降低生成图像的图像质量。

ct图像的模型主要由两部分组成，即有效的人体组织反射信号和无效的噪声信号，而噪声信号则包括乘性噪声和加性噪声，其中加性噪声相对于乘性噪声来说对ct图像的影响非常微小，在通常的处理中我们忽略加性噪声。这样一来ct电信号的通用模型可以表示为：

o(x,y)＝p(x,y)q(x,y)(1)

其中，x表示图像的横坐标，y表示图像的纵坐标，r(x,y)表示无噪声信号，n(x,y)表示相乘性噪声。

为了便于噪声的分离，需要对以上的模型进行对数变换将相乘的式(1)模型变换成相加模型：

log(o(x,y))＝log(p(x,y))+log(q(x,y))(2)

步骤2)对图像进行剪切波变换多尺度和多方向分解；

对上一步对数变换后便于噪声分离的ct图像应用多尺度分解，经分解后得到一张和原图等大的低频ct图像分量和多张同样与原图等大的高频ct图像分量。其中低频部分不做处理，在步骤3中会进行下一步处理，对于高频部分采用剪切滤波器作用到各个ct子带，对各尺度子带图像用剪切滤波器组进行方向分解。

步骤3)使用tv-l1模型分解图像为cartoon和texture部分；

运用tv-l1模型可以将图像分解成平滑的cartoon部分和保留纹理和噪声的texture部分，tv-l1模型可表示为：

其中u和f分别是输出图像和输入图像，tv-l1模型在不同的尺度下得到的结构不同，在尺度空间上能够有效的把不同的结构区分开。参数λ取值的大小对图像分解平滑有很大的作用，λ取值越小则平滑的越细致，cartoon部分的平滑效果非常干净。并且分解的另一部分，也就是细节纹理部分(texture)也非常清晰。这里取上一步中的低频部分和这里的cartoon部分相混合组成新的低频部分，为下一步的处理做好准备。图1是经典图像处理图片lena经tv-l1模型分解所得的cartoon部分和经分解处理所得的低频部分的混合图。

步骤4)对混合图应用快速非局部均值去噪；

对混合图应用利用积分图像技术加速得到的快速非局部均值去噪，得到新的低频子带。设含噪声图像为v，去噪后的图像为u。u中像素点x处的灰度大小的获取可以表示为：

其中权值w(x,y)表示像素点x和y之间的相似度，由以x，y为中心的矩形邻域v(x)、v(y)间的距离||v(x)-v(y)||²决定：

z(x)是归一化系数，h是平滑参数，h越大去噪程度越深，但如果过大则会导致图像变得模糊；如果h越小，越能够保留足够多的边缘细节成分，但h太小又会有过多噪声。

由于原始非局部均值算法复杂度太高，计算耗时非常大，有必要对其进行改进。

这里采用快速非局部均值方法主要由积分图像技术加速，积分图像技术是一种计算图像中矩形区域和的快速方法，积分图像技术的长处只要先计算出积分图像我们就能够在常量时间范围内计算图像任何大小的矩形的和。在积分图像(integralimage-ii)上任意一个位置(x,y)上的ii(x,y)表示该点左上角矩形范围内所有像素的和，表示如下：

其中i(x,y)表示输入图像上相关位置的像素点。若给定图像i，则从上到下、从左到右计算得到和的积分图像公式表示如下：

ii(x,y)＝ii(x-1,y)+i(x,y)+ii(x,y-1)-ii(x-1,y-1)(8)

其中(x＜0||y＜0)时ii(x,y)＝0,i(x,y)＝0。得到积分图像之后，图像中任意矩形区域和通过如下公式计算：

sum(m,n)＝ii(x,y)+ii(u,v)-ii(x,v)-ii(u,y)(9)

其中矩形大小为：m＝x-u,n＝y-v。

在利用积分图像技术为非局部均值去噪加速的应用中，假如待处理图像一共有像素点n个，设置搜索窗口为d×d(d＝2×ds+1)大小，设置邻域窗口为d×d(d＝2×ds+1)大小，那么两个矩形邻域之间的计算相似度时间就是o(d²)，其中每一个像素点都需要计算与搜索窗口内任何一个像素之间的相似度，而搜索窗口内具有d²个像素，所以可得非局部均值nl-means的计算复杂度就是o(nd²d²)。想要提高nl-means的运行计算速度就好从其计算复杂度o(nd²d²)入手，就该计算复杂度而言可以从o(d²)入手，也就是从降低邻域间相似度的计算入手。在原始计算过程中，想要计算邻域之间的距离的时候需要对每一个像素点求得差值，遍历所求距离的两个邻域。首先构造一个关于像素差值的积分图像可以解决这个问题：

其中st(x)＝||v(x)-v(x+t)||²，由此在计算两个邻域v(x)和v(y)(y＝x+t)之间的距离的时候，就可以在常量时间内完成计算：

为了提高运算速度，降低空间复杂度，以上所讲的算法把偏移量当做最外层的循环，也就是说每次只需要在一个偏移方向上求取积分图像，并对该积分图像进行处理即可。而不需要一次性求取出所有积分图像。由此，原计算复杂度o(nd²d²)降低为o(nd²)，加快了算法运行速度。

步骤5)对高频子带的剪切波系数进行阈值收缩处理；

在对高频部分进行阈值收缩处理的时候，其阈值函数的选取在很大程度上决定这去噪效果的好坏。本发明采用新型的自适应阈值收缩算法。对得到的多个高频ct子带的剪切波系数采取自适应阈值收缩处理，本发明使用的针对ct医学图像的阈值函数：

式中，σn是噪声的标准差，m是对应变换域内变换系数的总体个数，tj代表j层的自适应参数，经剪切波分解后不同分解层的系数会有不同的分布，所以tj基于不同的j层及具体实验目标选取；在本发明中，512*512的ct图分解为一个低频部分和4个高频部分，阈值t≈5σn，对于细节较多的分解子带，可使用t≈4σn或t≈3σn。

步骤6)对处理后的系数进行剪切波逆变换；

经过阈值收缩以及快速非局部均值处理有就得到了去噪后的剪切波系数，想要获得去噪之后的ct图像，就需要对处理后的系数进行剪切波逆变换，由此得到方便医师诊断分析的清晰的ct图像。

本发明整体步骤流程图如图2所示。

案例分析

本发明通过以具体的医学ct图像为对象，用基于快速非局部均值和tv-l1模型的剪切波变换医学ct图像去噪方法进行试验，高频子带中采用自适应的阈值算法，以psnr值为评价指标进行验证(psnr值越大代表去噪效果越高)，同时通过与现有技术对比展现本发明的优越性能。实验以医学ct噪声图像经典lena图(图片大小为512×512)作为输入数据带入不同去噪方法中进行试验，案例分析整体流程图如图3所示。实验通过对比fdct(快速离散曲波变换)，ffst(快速有限剪切变换)。各种算法应用在lena图的实验效果图以及原图和噪声图如图4，各种算法应用在头部ct图的实验效果图以及原图和噪声图如图5所示。

表1，2中可看出，从经典图像lena和医学ct图像的实验数据可以看出，噪声方差越大，对去噪算法的要求更高。本发明能够在噪声方差增大的情况下，保持相对稳定的去噪效果。在同一噪声方差上，效果也均好于fdct和ffst。

表1：lena图不同去噪算法在不同噪声的psnr/db值

表2：医学ct图不同去噪算法在不同噪声的psnr/db值

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的同等技术手段。