医学图像变换

1、第四章 医学图像变换在医学图像处理与分析中广泛应用着各种图像变换技术，它们是图像处理与分析的重要工具之一，通过各种图像变换来转换图像的表示域以及表示数据，给后续的图像处理工作带来极大的方便。图像变换是一种为了达到某种目的而对图像使用的一种数学操作，经过图像变换后的图像将能够更方便、更容易地被处理和操作，因此图像变换在图像增强、图像复原、图像编码、特征抽取等方面有着广泛的应用。例如，傅立叶变换可使处理分析在频域中进行，使运算更简便；某些图像经过变换后往往能反映出图像的灰度结构特征，从而更便于分析；还有许多变换可使变换后的能量集中在少数数据上，从而便于实现数据压缩、图像传输和存储等等方面。在实际的

2、图像处理中，图像变换可以看作是一个数学问题，即对原图像函数寻找一个合适的变换核，但本质上来说，图像变换有着深刻的物理背景。常用的图像变换方法主要有：傅立叶变换、余弦变换、小波变换、哈达玛变换、KL变换、哈尔变换、斜变换等。由于傅立叶变换和小波变换目前应用的较为普遍，并且在理论上也比较重要，所以本章将重点讨论这两种图像变换形式。第一节 傅立叶变换傅立叶变换是一种正交变换，它广泛地应用于很多领域，从某种意义上说，傅立叶变换就是函数的第二种描述语言，掌握了傅立叶变换，人们就可以在空域和频域中同时思考处理问题的方法。由于它不仅能把空间域中复杂的卷积运算转化为频率域中的乘积运算，还能在频率域中简单而有效

3、地实现增强处理和进行特征抽取，因而在图像处理中也得到了广泛的应用。 一、一维傅立叶变换一维连续信号的傅立叶正变换和反变换的数学表达式如下： （4.1） （4.2）从上式可以看出F(u)通常是自变量u的复函数，可以将其表达为如下形式： （4.3）可以得到： ， （4.4） （4.5）其中称为的频幅谱或傅立叶谱，为傅立叶变换的相位角，振幅谱的平方通常被称为的能量谱。傅立叶变换的物理含义是将信号f(x)表达为一系列正交基函数的加权求和。傅立叶正变换的目的就是求出对应各正交基的权重。傅立叶反变换的目的就是通过这些基函数的加权和，恢复出原始信号。同样道理，在离散域中，一维连续傅立叶变换式中的积分运算可以

4、简化为求和，一维离散傅立叶正变换和反变换的表达式如下： （4.6） （4.7）如下图所示是一个淹没在噪声当中的信号，仅从时域信号看，很难看到信号本身的特征。将该信号进行傅立叶变换后并如图显示所示，则能很容易看出信号的主要信息特征。图4-1 信号的时域显示图4-2 信号的频域显示二、二维傅立叶变换对于二维信号或者图像信号来说，同样存在连续傅立叶变换和离散傅立叶变换，它们的正变换和反变换分别如下式所表示：二维连续傅立叶变换： （4.8） （4.9）二维离散傅立叶变换： （4.10） （4.11）三、傅立叶变换的性质（1）平均值傅立叶变换域原点的频谱分量是空间域的平均值的倍，即式（4.12）所示：

5、（4.12）（2）变换的周期性设m，n为整数，m，n=0，1，2，将u+mN和v+nN代入式中右边，有： （4.13）上式中，右边第二个指数项为单位值，因此傅立叶变换是周期性的，即：（3）对称共轭性由离散傅立叶变换定义可方便地证明，傅立叶变换满足： （4.14）（4）平移性如果用表示傅立叶变换对，则平移性是指： （4.15） （4.16）由于，因而说明的移动并不影响它的傅立叶变换的幅度。（5）线性特性和比例性设，若和分别是的傅氏变换，则根据定义可知其线性特性为： （4.17）同时，也容易证明其比例性为： （4.18）（6）可分离性可分离性的主要优点是可以通过两次一维变换来实现一个二维变换（或反

6、变换）。可分为两步，第一步：先沿轴对求一维离散傅立叶变换，得中间结果，即： （4.19）第二步，再沿轴对求一维离散傅氏变换，得最后结果，即： （4.20）（7）旋转性质由于在极坐标下表示二维函数图形的旋转特性非常方便，所以可以将坐标进行转换。空间域坐标变换为： ， （4.21）频率域坐标变换为： ， （4.22）便是极坐标中的傅氏变换对。可以证明二维离散傅氏变换具有如下旋转性质： （4.23）（8）微分性质傅立叶变换的微分性质可表示为： 和 （4.24）作为特例，在图像增强中用到的拉普拉斯（Laplace）算式，可以定义为： （4.25）则由微分性质可知laplace算子的傅氏变换为，即：便是

7、在模式识别技术中经常用到的laplace算子。四、快速傅立叶变换（FFT）对于一维的傅立叶变换的计算来说，每次计算都需要进行N次复数乘法和N1次加法，用这个公式计算N点的一维傅立叶变换就有计算次的乘法和加法运算。而快速傅立叶变换把只计算一次，然后把它存放在一个表里，以备查用，这样可以使计算次数减少为，因此原始变换算法与快速傅立叶变换算法的计算量之比为，当N比较大时，计算量节省是相当可观的。下面介绍一种逐次加倍法的快速傅立叶变换算法：令： 其中N为2的正整数次幂，且，则上式可表示为： （4.26）由于性质，所以可以将上式写为： （4.27）现在可以定义： （4.28） （4.29）因此可以将（4

8、.1.26）简化为： （4.30）同理，由和可得： （4.31）现在我们来仔细分析式（4.28）到式（4.29），式（4.30）和式（4.31）表明1个N点的变换可通过将原始表达式分成两半来计算，对第一半的计算需要根据式式（4.1.27）和式（4.1.28）计算2个（N/2）点的变换，这样所得到的所得到的可以代入式（4.30）和式（4.31）以得到时的。对剩余的的计算也与此类似。因此该算法的关键是将输入数据排列成满足连续运用式（4.28）和式（4.29）的次序。例如要计算1个8点的快速傅立叶变换，需要将它们排列成的形式，如下图所示，这样在第1层先计算4个2点变换，在第2层用以上4个结果计算2个

9、4点变换，在第3层再用医师2个结果计算1个8点变换。图4-3 快速傅立叶变换的输入数据排序对输入数据的排序可根据一个简单的位对换规则进行，如用表示的1个自变量值，那么它排序后对应的值可通过把表示成二进制，并左右对换各位得到。例如时，排序后为，因为，而反过来的。因此，只有把输入数据进行重新排序，则输出结果才是正确的次序，对于傅立叶反变换也可以类似计算。五、傅立叶变换的物理意义傅立叶变换的本质是将一个时域上的信号转换到频率域。我们平时观测到的各种物理量都是随时间的变化，进行傅立叶变换以后就是这些物理量在不同频率上的分布。从傅立叶变换原理可以看出，任何连续测量的时序信号，都可以表示为不同频率的正弦波

10、信号的无限叠加。对于图像傅立叶变换的物理意义来说，是将以灰度信息表示的图像转变成以不同频率信息表示的图像，换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，而图像的频率信息是表征图像中灰度变化剧烈程度的指标。因此，和对信号的傅立叶变换相似，任何图像的傅立叶变换也可以用许多基图像的线性组合来表示，如图4.4，图4.5，图4.6所示。uv图4-4 基函数图像图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率

11、值较高。傅立叶变换以前，图像是由对在连续空间上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度？因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。从傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小。一般来讲，梯度大则该点的亮度强，否

12、则该点亮度弱。这样通过观察傅立叶变换后的频谱图，如图4.5和图4.6所示，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。图4-5 以数字矩阵形式表示的图像傅立叶变换图4-6 以基函数表示图像傅立叶变换六、基于傅立叶变换的频域滤波空间域和频率域的线性滤波基础都是卷积定理，该定理可以表示为: （4.32） （4.33）其中，符号*表示两个函数的卷积，双箭头两边的表达式组成了傅立叶变换对。第一个表达式表明两个空间函数的卷积可以通

13、过计算两个傅立叶变换函数的乘积的逆变换得到。相反，两个空间函数的卷积的傅立叶变换恰好等于两个函数的傅立叶变换的乘积。同样的情况也出现在第二个表达式中。基于傅立叶变换的频域滤波问题，主要通过第一个表达式来实现。在空间域中的操作图像与设定滤波函数。根据卷积定理，我们可以在频域中通过乘以来得到相同的结果，我们通常将空间滤波器的傅立叶变换称为滤波传递函数。 基本上，频域滤波的目的是选择一个滤波器传递函数，以便按照指定的方式修改。基于傅立叶变换的频域滤波处理的过程，包括如下主要环节：图4-7 频域滤波基本过程根据以上分析，针对一幅医学图像进行频域滤波处理，其结果如图4.8和图4.9所示：图4-8 医学图

14、像及其频谱显示图4-9 滤波器图形及其滤波结果七、 傅立叶变换相关函数在MATLAB中较为常见的傅立叶变换相关函数主要有以下这些：傅立叶变换相关函数conv卷积conv22维卷积fft快速傅立叶变换fft22维快速傅立叶变换ifft逆快速傅立叶变换ifft22维逆快速傅立叶变换filter离散时间滤波器filter22维离散时间滤波器abs幅值angle四个象限的相角unwrap在360边界清除相角突变fftshift把FFT结果平移到负频率上表4-1 傅立叶变换相关函数第二节 离散余弦变换离散余弦变换（DCT）是利用傅立叶变换的对称性，采用图像边界褶翻操作将图像变换为偶函数形式，然后对这样的

15、图像进行二维离散傅立叶变换，变换后的结果将仅包含余弦项，所以称为离散余弦变换。离散余弦变换可以将图像描述为不同幅值和频率的正弦值之和的形式。离散余弦变换的一个重要的性质就是，有关图像的许多重要可视信息都集中在变换后的一小部分系数当中，因此，在图像压缩处理中离散余弦变换是非常有用的，而且离散余弦变换已经是有损压缩国际标准JPEG算法的核心变换。一、 一维离散余弦变换如果首先定义一个常量为： （4.34）那么一维离散余弦变换（DCT）的正变换核由下述关系式给出： （4.35）同样一维离散余弦的反变换也可以定义为： （4.36）二、二维离散余弦变换同样道理，二维离散余弦变换的正反变换相同也可以表示为

16、： （4.37） （4.38）不难看出，二维离散余弦变换的变换核是可分离的，因而可通过两次一维变换实现一个二维变换，二维余弦变换在图像压缩编码中得到了广泛应用。如图4.10所示，图像经DCT后，能量集中于频率平面的左上角，因此DCT常常用于医学图像数据的压缩。DCT图4-10 离散余弦变换频谱显示三、离散余弦变换的MATLAB实现在MATLAB的图像处理工具箱内，提供了两种不同的离散余弦变换的计算方法。第一种方法是直接使用dct2函数，该函数是一个基于FFT算法来提高当输入矩阵较大时的计算速度，dct2函数的使用格式如下：B=dct2( A，M，N )；其中A表示要变换的图像，M、N是可选参数

17、，表示填充后的图像大小，B表示变换得到的图像矩阵，具体程序如下所示：A=imread(mri2.jpg);B = dct2(A);subplot(1,2,1), imshow(A);subplot(1,2,2), imshow( log(abs(B); 图4-11 离散余弦变换前后的图像显示第二种方法是使用函数dctmtx得到离散余弦变换的矩阵，这种方法适合于较小的输入方阵，其使用格式为B=dctmtx（A），其中，A表示要变换的矩阵，B表示变换后得到的矩阵。四、离散余弦变换的初步应用离散余弦变换是JPEG压缩算法的基础，下面对离散余弦变换在图像压缩中的具体应用做一个简单的介绍，在JPEG图像

18、压缩算法中，输入图像被分为8*8或者16*16块，然后对其中每一个块进行二维离散余弦变换，接着再对离散余弦变换的系数进行量化、编码和发送，JPEG的接收者只要通过对离散余弦变换系数进行解码，计算每个块的二维逆变换，然后将所有的块重新组装成一幅图像。对于一般的图像来说，离散余弦变换系数当中许多都是接近于0的数值，可以丢弃这些系数而不会对图像的重建质量产生重大影响，下面就以对一幅医学图像的压缩处理为例来作为说明。A=imread(mri2.jpg);B = dct2(A);B(abs(B) 0，a=1 （c）b不变，a=2图4-14 小波伸缩示意图小波变换之所以在时频域中分辨率都可随频率变化而变化

19、，就因为在高频时使用小尺度a值，时轴上观察范围小，而频域上相当于用高频小波作细致观察；低频时使用大尺度a，时轴上考察范围大，而频域上相当于用低频小波作概貌观察。虽然分析频率有高有低，但在各分析频率段内分析的品质因数却保持一致，这是一项很符合实际工作需要的特点，因为如果希望在时域上观察得愈细致，就愈要压缩观察范围，并提高分析频率。利用小波变换所具有的这种数学显微镜特点和频域带通特性，可以把所需的信号分离出来，进行分析研究。（4）平移（Shifting）对小波作平移就是简单地对原信号进行延时。在数学上，对函数f(t)延时k可表示为f(t-k)：图4-15 小波平移示意图从小波函数和伸缩函数的一些特

20、性上，我们可以来区分二者的关系，在整个实空间内，对的积分为零，即，因此用来定义“细节”部分，但对的积分为1，即，因此用来定义“近似”部分，小波子空间横跨相邻的二个伸缩子空间：，各小波子空间之间为正交的，不互相包含，即。二、 连续小波变换：连续小波变换可定义为： （4.42）这里从上式可知，小波变换是信号与小波基函数的内积，是对信号满足一定附加条件的滤波，这种附加条件反映在小波函数及伸缩因子的选择上。由于变量s，t是连续变量，所以把这种小波变换称为连续小波变换（continuous wavelet transform，CWT）对连续小波进行逆变换可获得原信号，小波逆变换的表达式为： （4.43）

21、这里：连续小波变换具有以下重要性质：(1) 线性：一个多分量信号的小波变换等于各个分量的小波变换之和；(2) 平移不变性：若的小波变换为，则的小波变换为；(3) 伸缩共变性：若的小波变换为，则的小波变换为；(4) 自相似性：对应不同尺度参数s和不同平移参数b的连续小波变换之间是自相似的；(5) 冗余性：连续小波变换中存在信息表述的冗余度。三、 离散小波变换离散小波变换可定义为： （4.44） （4.45）并且其逆变换为： （4.46）其中，和分别称为“近似系数”和“细节系数”。信号的离散小波变换过程如下图所示，将信号分解成近似系数与细节系数：图4-16 信号小波分解示意图这里，如果滤波器的长度

22、为2N，信号的长度为n，则信号F和G的长度为n+2N-1，则和的长度为： 同样道理，我们可进一步将近似系数按上面同样的方式分解为和，并可以依次循环下去（如下图所示），直到得到满足需要的结果为止。图4-17 近似系数的再分解当分解层次J=3时，可用树形结构来清晰表示这个分解过程：图4-18 三层小波分解示意图而离散小波的逆变换，就相当于是对信号的合成过程，合成过程是信号分解的逆过程，通过插零值（上采样）及对插值结果滤波，来重建上一级的近似细节信号，并最终可以得到合成信号。一维离散小波逆变换的过程如下图所示：图4-19 一维离散小波逆变换四、快速小波变换像快速傅立叶变换FFT一样，快速小波变换利用

23、小波系数的自相似特性来实现离散小波变换的快速计算，它也被称为Mallat算法，即：某一分辨率下，离散小波变换的“细节系数”和“近似系数”，都可以通过是上一级高分辨空间中的“近似系数”和小波向量的卷积，而后对卷积结果进行亚取样，从而计算得到。它们之间的运算关系可以用下列关系式表达： （4.47） （4.48）其中，卷积只在非负偶数时进行计算，所以与以2为步长进行抽样的效果相同。五、 二维小波变换从一维小波变换，很容易推广到二维小波变换，这里主要包含一个2D的伸缩函数（表示近似部分），和三个2D的小波函数。 伸缩函数 测量水平方向上的细节变化 测量垂直方向上的细节变化 测量对角线方向上的细节变化对

24、j级的近似部分进行2D小波变换，可得到下一分辨级（j+1级）的一个近似和三个方向的细节（水平、垂直、对角线），对这个j+1级近似仍可继续不断地分解。二维小波变换主要应用于对二维图像的处理中，用2D小波变换对图像的分解和重建过程如下图所示：图4-20 二维离散小波变换图4-21 二维离散小波逆变换分解层次J=2时，用树结构的形式表示成： 图4-22 二维小波分解结构树形示意图 图4-23 一幅磁共振图像的二层小波分解效果六、 常用小波函数及其特性一般要根据实际问题的需要，选择和构造不同的小波，但实际上，小波函数很难建立，因此我们往往直接利用现有的小波函数。下面介绍几个常用的著名小波家族及其性质。

25、小波函数HaarDaubechiesBiorthgonalCoifletsSymletsMorletMexican hatMeyer小波缩写名baardbbiorcoifsymmorlmexhmeyr表示形式haardbNbiorNr.NdcoifNsymNmorlmexhmeyr举例baardb3bior2.4coif3sym2morlmexhmeyr正交性有有无有有无无有双正交性有有有有有无无无紧支撑性有有有有有无无无连续变换可以可以可以可以可以可以可以可以离散变换可以可以可以可以可以不可以不可以可以支撑长度12N-12Nr+12Nd+16N-12N-1有限有限有限滤波器长度22Nmax(

26、2Nr,2Nd)+26N2N-4,4-5,5-8,8对称性对称近似对称不对称近似对称近似对称对称对称对称表4-2 常用小波函数及其特性七、小波变换与医学图像处理数字图像处理的发展开始于20世纪60年代初期，至今已有40多年研究历史，其经典的图像处理方法有很多。而小波变换则被看作是一种用于多层次图像分解的数学工具。图像数据经过小波变换后可以用小波系数来描述，小波系数体现出原图像数据的性质，图像数据的局部特征可以通过对小波系数进行处理而改变。小波变换在医学图像处理上的应用思路主要采用将空间域上的图像数据变换到小波域上，得到多层次的小波系数，根据所采用的小波基特性，分析小波系数特点，针对不同需求，即

27、可以结合常规的图像处理方法，也可以提出更符合小波分析的新方法来处理小波系数，然后再对处理后的小波系数进行小波逆变换，得到所需的目标图像。因此，小波分析与变换在医学图像恢复、图像增强、图像分割、图像数据库检索、图像配准和融合、图像重建等处理上都能得到应用。基于小波变换的医学图像处理过程可以用下图表示。小波正变换图像处理小波逆变换图像输入图像输出图4-24 小波变换与图像处理（1）基于小波变换的阈值化去噪方法传统的去噪方法是将被噪声干扰的信号通过一个滤波器，直接滤掉噪声频率成分，但对于脉冲信号、白噪声、非平稳过程信号等，传统方法存在一定局限性，对这类信号，在低信噪比情况下，经过滤波器处理，不仅信噪

28、比得不到较大改善，而且信号的位置信息也被模糊掉了。基于小波变换的去噪方法利用小波变换中的变尺度特性对确定信号具有一种“集中”的能力。如果一个信号的能量集中于小波变换域少数系数上，那么对这些系数的取值必然大于在小波变换域内能量分散于大量小波系数上的信号或噪声的小波系数值。因此，基于小波变换的域值化去噪方法，关键是对小波变换系数如何进行筛选，它依赖于对小波系数的域值化和域值门限的选取。小波域值化去噪方法就是在小波系数进行取舍之前，实际上按着一定准则，将小波系数划分成两类：一类是重要的、规则的小波系数；另一类被看作是非重要的或者受噪声干扰较大的小波系数。通常以小波系数的绝对值作为小波系数的分类单元，

29、小波系数绝对值趋向零，意味着小波系数所包含的信息量少，并且强烈地受噪声干扰。因此给定一个小波系数阈值，所有绝对值小于某个阈值的小波系数被划为“噪声”，它们的数值用零代替；而超过阈值的小波系数的数值用阈值缩减后再重新取值，这种方法意味着阈值化或者缩减小波变换将在小波域中移去小幅度的噪声或者非期望的信号。最后再进行小波逆变换就得到所需要的恢复图像了。 （b）含噪声的原始图像 （b）用小波变换方法去除噪声的结果图4-25 小波变换与图像去噪（2）小波变换与图像增强图像增强技术的目的是采用一些技术手段，有选择地突出图像中感兴趣的特征或者抑制图像中不需要的特征，改善后的图像比原始图像更加满足某些特殊分析

30、的需要。图像增强处理主要有两类方法：空域法和频域法，而频域增强的关键就是要借助具体的频域变换，来实现图像信息的转换，频域变换一般指傅立叶变换，但它也可以延伸为其他变换，如DCT变换、Walsh变换和小波变换。在图像增强中，可以用来反映图像增强前后之间的函数关系，这种函数关系可以是线性的，也可以是非线性的。虽然小波变换在图像增强的应用中仅充当一个频域变换的作用，但是增强算法的设计可以充分利用小波分析的时一频局部化特性，更加有效地提高图像增强的质量和算法的时效性。比较常用的小波增强方法有：子带增强法、反锐化掩模法、自适应增益增强法等等。 （a）原始核磁共振图像 （b）小波变换增强结果图4-26 小

31、波变换图像增强（3）小波变换与图像编码最新的小波编码方法大多是沿袭以前变换编码的基本思想，即去相关性，因此变换、量化和熵编码是构成小波变换图像编码的三个主要成分，目前，人们主要围绕这三个方面的内容开展小波变换图像编码的研究。小波编码在实现结构和方法上与子带编码非常一致，可以说小波编码是子带编码的特例。一般来说，子带编码并非强求使用完全重构条件，但用小波编码时则使用严格的完全重构滤波器，并且要求符合正则性条件。小波编码的基本思想：将原始图像经二维小波变换后，转换成小波域上的小波系数，然后对小波系数进行量化编码。由于小波变换后使得原始图像能量集中在少数部分的小波系数上，因此最简单的系数量化方法就是

32、将某一阈值以下的系数略去，或者表示为恒定常数，只保留那些能量较大的小波系数，从而达到数据压缩的目的。在这里，可以采用标量量化或矢量量化的方法分别在不同分辨率的小波频带上来完成的。 （a）原始核磁共振图像 （b）小波编码压缩结果图4.27 小波变换与图像压缩编码上图对一幅磁共振图像进行了小波编码压缩，采用sym4进行3层小波分解，选择域值为24.18，采取全局域值压缩方法，压缩结果为能量保留99.47，小波系数置零率为91.63，图b是经过解压缩后的图像，可以看出整个图像质量和原始图像相差不大。（4）小波变换与图像融合医学图像融合的方法大致可以分为两类，一是直接在像素域上的融合，二是基于变换域的融合。在医学图像融合中,目前所用的像素域方法有加权叠加法、色度空间融合法等，变