基于多尺度Retinex的医学CT图像增强算法

曹 莉，许玉龙

(河南中医药大学信息技术学院，河南 郑州 450000)

1 引言

图像增强技术是通过某种手段变换或增加原图像数据，从而改善图像视觉效果，使图像匹配视觉特性[1]。医学CT图像一般表现为信息量大且图像噪声多而杂、图像暗区模糊等，严重影响医学CT图像的病理分析，从而导致医生临床诊断结果产生较大的偏差，因此对医学CT图像进行图像增强就变得尤为重要。

目前大量学者对医学图像增强方法进行了研究，并取得了一定的研究成果。文献[2]提出基于Canny算子加权引导滤波的Retinex医学图像增强算法。采用Canny边缘检测算子和加权引导滤波器，估计边缘权重和光照，获取入射光分量，根据Retinex算法，计算反射光分量，量化处理图像后输出增强图像，该算法能够保留图像细节信息，但是图像灰度分布范围较小。文献[3]提出基于剪切波和改进Pal-King的图像增强算法。采用剪切波变换分解图像，通过自适应阈值去噪去除图像噪声，利用剪切波反变换重构图像，运用Pal-King算法增强图像。该算法能够有效改善图像的视觉呈现效果，但是该方法的图像边缘信息不完整。

针对上述问题，提出基于多尺度Retinex的医学CT图像增强算法。通过提取医学CT图像暗区信息，获取医学CT暗区图像，线性平滑处理医学CT暗区图像，采用多尺度Retinex算法，计算医学CT图像像素间的明暗关系，转换后得到增强医学CT图像。所提算法能够扩大灰度分布范围的同时，完整保留图像边缘信息，医学CT图像增强视觉效果较好。

2 医学CT暗区图像获取

恢复医学CT图像的关键就是获得准确的光照分布，首先将图像转换到RGB彩色空间[4]，获得图像的RGB三通道图像，然后将医学CT图像描述为

I(x)=J(x)t(x)+A(1-t(x))

(1)

式(1)中，I(x)表示医学CT图像；J(x)表示期望获得的辐射值；A表示光照强度分布；t(x)表示没被散射的光线部分。提取医学CT图像暗区信息，定义医学CT图像I(x)中每个像素点的三通道最小值构成的图像为Imin

(2)

对应的图像B为

(3)

式(3)中，Q(i，j)是以像素Ic(i，j)为中心的图像块，(x，y)表示Q(i，j)内的任意像素点。

若医学CT图像I(x)为室内光照下图像，将B称为医学CT图像I(x)的暗区图像[5]，室内环境下会存在很多阴影或彩色，导致医学CT图像的一些信息量较低，因此以此为理论依据可估计光照强度分布如下

B(i，j)=Jdark(i，j)t(x)+A(i，j)(1-t(x))

(4)

式(4)中，Jdark(i，j)表示平均照度，因Jdark(i，j)的值无限接近于零可以忽略，则式(4)等价于

B(i，j)=A(i，j)(1-t(x))

(5)

故可以获得医学CT暗区图像的光照强度部分图A(i，j)可表示为

(6)

综上可得在图像中，光照强度分布图与暗区图像的大小和系数1-t(x)相关，由此获取医学CT暗区图像。

3 医学CT图像线性平滑处理

由于获取的医学CT暗区图像是通过计算以其自身为中心的窗口内最小值获得的，但是若某一像素点的某一通道值极小时，所对应的像素点都将成为这一极小值，就产生了分块特性[6]。而光照强度通常是平滑过渡的，因此为解决这个问题还需要对其进行线性平滑处理[7]。

线性平滑处理还需要满足以下两个要求：

1)医学CT暗区图像中灰度值没有产生明显变化的低频区域；

2)保持医学CT暗区图像中灰度值变化较为明显的高频区域，以确保能够完整地保留图像边缘信息。

为完成以上的线性平滑处理要求，在以医学CT图像I(x)上的一点为中心，大小为M×M的窗口Wk(k=1，2，…，w×h)内的每个像素通过式(7)进行线性近似

E=akB+bk

(7)

式(7)中，ak、bk为近似参数。

M的大小由图像的大小w×h来决定，E的结果会随着近似参数的变化而变化，E和B之间的差可由式(8)求得

(8)

式(8)中，ε为数值较低的常数用来平衡ak防止其过大。

使用线性回归对R(ak，bk)进行求解如下

(9)

式(9)中，uk和δk分别为B的平均值和方差。Ccount(W)表示窗口内像素数量。

当R(ak，bk)取最小值时，在低频区域内，B和uk的值近似相等。线性近似参数的求解结果近似于ak=0，bk=uk，E=uk，即对窗口内全部像素的灰度值平均化；在高频区域内ak、bk的解近似于ak=1，bk=0，E=B，即保持窗口内原有灰度值[8]以满足医学CT暗区图像线性近似的要求。

求得所有ak，bk以后，由于可能存在一个像素会同时在多个窗内出现的问题，所以在获得平滑结果时，取窗Wk内全部ak、bk的平均值进行运算，即：

在对医学CT暗区图像进行线性平滑处理后，不但解决了分块效应的问题，并且修正了出现错误的区域，获得用于医学CT暗区图像增强的光照强度分布图。

4 基于多尺度Retinex的医学CT图像增强

多尺度Retinex算法认为人眼处理输入信息是分别按照长波、中波、短波单独进行处理的。由于部分医学CT图像为彩色图像，因此也需要将其分解为长波、中波、短波三个波段与图像R，G，B中三个分量相互对应。为获得分解后的医学CT图像将像素中R分量归为第一幅图像，G分量归为第二幅图像，B分量归为第三幅图像，三幅图像分别对应RGB三个波段的反射光强度。由于处理像素间明暗关系会使用到浮点运算。为保证后续的运算精度，需要将整数格式的像素值转换为浮点数。

单尺度Retinex算法[9]的计算模型可以表示为

R(x，y)=logS(x，y)-log(F(x，y)⊗S(x，y))

(10)

式(10)中，R(x，y)代表增强医学CT图像，S(x，y)代表未被处理的医学CT图像，F(x，y)为环绕函数，其表达式为

F(x，y)=λ*exp(-(x2+y2)/2c2)

(11)

式(11)中，λ表示尺度参数，c表示环绕尺度，其数值大小对医学CT图像处理质量有直接影响。但是单尺度Retinex算法处理后医学CT图像清晰度较低，并且会出现光晕，因此对加权求和与单尺度Retinex模型进行融合，得到多尺度Retinex算法模型，如图1所示。

图1 多尺度Retinex算法模型示意图

多尺度Retinex算法模型可以表示为

RMSRI(x，y)

=Qk{logSi(x，y)-log(Fk(x，y)⊗Si(x，y))}

(12)

式(12)中，RMSRI(x，y)为增强后的Si(x，y)通道输出，k表示尺度数量，Qk表示k的相关权值，Fk(x，y)为高斯环绕函数[10]。利用式(12)来提高图片色彩恢复效果

Ci(x，y)=β[αlogSi(x，y)]

(13)

式(13)中，β为增益常数，α为线性强度。通过对参数的反复研究得出当α=123，β=45时，医学CT图像的最终处理效果最好，保留原图中通道之间的相互关系，使最终获得医学CT图像的色彩和实景场景更为接近。通过修正获得的输出显示为

rMSRI(x，y)

=Qk{logSi(x，y)-log(Ci(x，y)⊗Si(x，y))}

(14)

式(14)中，rMSRI(x，y)表示修正获得的输出。由于经多尺度Retinex算法处理后，原医学CT图像中最亮和最暗部分都会不同程度接近灰色，导致增强后的医学CT图像亮度范围变小，真实感变差[11]。因此对原多尺度Retinex算法进行如下处理

r(x，y)=log(S(x，y))-λlog(L(x，y))

(15)

式(15)中，L(x，y)表示归一化光照估计值。

多尺度Retinex算法是基于假设为先验来估计光照分量的，因此求得的反射分量并非是理想的增强医学CT图像，需对其明暗关系进行进一步处理。

若从医学CT图像中的(x1，y1)到(x2，y2)只有一条路径，该路径上存在N个像素，其值分别为(d1，d2，…，dn)。用连续相乘计算该路径起点S(x1，y1)和终点D(x2，y2)之间的明暗关系，可表示为

(16)

式(16)中，T(x)为门限函数，也就是说若相邻像素值变化低于既定门限时，即可认定其值没有发生变化。

由于式(16)需要取连续乘积的对数，因此首先可以取各个比值的对数，然后将所有结果进行加和，即

logD=T(x)(logd3-logd2)+…

+T(x)(logdn-logdn-1)

(17)

从式(17)可以看出，若将像素值d的对数logd直接相减也相当于计算其比值，即logd1-logd2=logd1/d2，同理连续相乘也可以通过将每个比值取对数相加来获得。因此得到可以通过连续相减或者连续相加的方法来计算医学CT图像像素间的明暗关系。

在像素明暗关系处理的过程中还存在两个问题需要解决，一是门限的判决问题，即当相邻像素的亮度变化低于既定门限时，就可以将其反射能力判定为没有发生变化[12]；二是起点与终点之间存在多路径问题，即在起点相同的情况下，像素点可通过不同的运动路径到达终点，选择路径不同所得到的结果也会随之发生不同。因此为解决以上问题，对医学CT图像进行如下处理：

设医学CT图像大小为m×n，并且初始医学CT图像的各点明暗值相同，即

lightness[i][j]=constant≤i≤m，1≤j≤n

(18)

首先对水平方向距离h=m/2的任意两个像素间明暗关系进行计算如下

(19)

然后计算垂直方向两个距离v=n/2的像素间明暗关系为

(20)

明暗值计算完成后，将水平方向距离调成h=m/4，垂直方向距离调成v=n/4，重复以上计算，直到两个方向的距离均为1时，完成一次迭代。利用上述算法分别对长、中、短波进行处理，获得三个波段内之间的明暗关系。根据多尺度Retinex算法，医学CT图像某点的色彩取决于三个波段的明暗值。若每个波段内的明暗为8比特，最亮部分为255，最暗部分为0，采用下式进行线性映射

(21)

式(22)中，min为在该波段内明暗值最小值，max为在该波段内明暗值最大值。Lb表示长波明暗值映射后的结果，转换后得到线性增强医学CT图像，由此实现医学CT图像增强。

5 仿真研究

为了基于多尺度Retinex的医学CT图像增强算法的有效性，在MATLAB6.1的软件环境下，将所提算法获得的医学CT图像灰度值和增强视觉效果与文献[2]算法、文献[3]算法进行对比分析。医学CT图像大小为512×512，使用OpenCV计算机视觉库的基本数据结构，以及一些图像操作的基本函数如cvLoadlmage、cvClonelmage等对图像文件进行打开及复制操作，以简化过程。

为了验证所提算法的医学CT图像增强效果，分别采用文献[2]算法、文献[3]算法与所提算法增强后的医学CT图像增强结果如图2所示。

图2 不同方法的医学CT图像增强效果对比结果

从图2中可以看出，采用文献[2]算法增强后的医学CT图像较为清晰，但仍存在失真现象，医学CT图像不够平滑；文献[3]算法由于在一定程度上合并了图像灰度级，导致医学CT图像较为模糊，缺少边缘信息；而所提算法增强后的医学CT图像较为平滑，视觉效果较好，能够有效保留图像边缘轮廓信息，并且明显改善图像失真的问题。因为所提算法线性平滑处理医学CT暗区图像，采用多尺度Retinex算法，计算医学CT图像像素间的明暗关系，使最终获得医学CT图像的色彩和实景场景更为接近，保持图像灰度值变化较为明显的高频区域，以确保能够完整地保留图像边缘信息，从而增强了医学CT图像视觉效果。

为了进一步验证所提算法的医学CT图像增强效果，分别采用文献[2]算法、文献[3]算法与所提算法的获取医学CT图像灰度值结果如图3所示。

图3 不同方法的医学CT图像灰度值对比结果

从图3可以看出，文献[2]算法和文献[3]算法增强后的医学CT图像灰度分布范围较小，而所提算法增强后的医学CT图像灰度分布范围较大。因为所提算法考虑分块特性，线性平滑处理医学CT图像，使医学CT暗区图像中灰度值未产生明显变化的低频区域，保持窗口内原有灰度值，以满足医学CT暗区图像线性近似的要求，从而提高暗区图像的信息量和对比度，扩大医学CT图像灰度分布范围。

6 结论

为实现医学CT暗区图像增强，扩大灰度分布范围，提高医学CT图像的视觉效果。本文通过提取医学CT图像暗区信息，获取医学CT暗区图像进行线性平滑处理，采用多尺度Retinex算法，获取图像中像素间明暗关系，对像素空间进行线性映射，得到增强医学CT图像。所提算法能够有效还原真实场景信息，扩大灰度分布范围，有效保留医学CT图像边缘信息。但是由于运算量较大，导致运行速度较慢，还需要在实际应用中对各个数据进行进一步优化。