基于梯度场变换的MRI图像增强方法

吴梦飞， 薛旭成， 兰太吉， 徐鑫伟

1. 中国科学院长春光学精密机械与物理研究所， 吉林 长春 130033；2. 中国科学院大学， 北京 100049

1 引言

MRI的全称是磁共振成像(magnetic resonance imaging)，在日常医学诊断中，MRI常用于肿瘤、炎症、退行性病变等疾病的检查。然而，由于磁共振机器在工作时会产生非常强大的磁场，其所成像的MRI图像存在较多的系统噪声，而且部分MRI图像整体昏暗，细节特征表现不清。这无疑给后续的诊断工作带来了很大的困难。针对于这一问题，目前有学者提出各种方法用于MRI增强，但这些方法在增强的同时也会对系统噪声产生过度增强的效果[1-4]。本文提出一种基于梯度场变换的图像增强方法，该方法在对MRI增强的同时可以有效地避免图像背景噪声的增强。实验结果表明，本文提出的方法增强效果优于直方图均衡化等传统的方法。

直方图均衡化是一种常用的图像增强算法，它能够将图像区域中的灰度范围进行拉伸，使其以相同的概率分布在整个直方图上，在医学影像领域得到了广泛应用。对于对比度较低的低照度图像，直方图均衡化技术能够有效地提升其对比度，产生较好的视觉效果。另外，因为直方图均衡化在理论上等同于熵最大化，它是熵最大化理论的直接体现，这意味着图像中的特征信息会以相同的概率表现出来[5]。但直方图均衡化也有其局限性，即容易对图像产生过度增强的效应。经过直方图均衡化处理过的图像，其部分细节特征会产生过亮的现象，同时在图像中灰度数量占比较少的微弱信息也会因为灰度级的拉伸而被合并[6]。图像伽马变换是另一种常用的空域图像增强方法。伽马校正可以有效地提升图像的对比度，但是对于低照度图像，普通的伽马校正难以产生较好的效果。文献[7]中提出了一种基于全局亮度平均值信息来对图像局部区域进行修正的伽马校正图像增强算法，该算法将图像分为高通滤波和低通滤波两个部分分别进行操作，通过融合图像的全局信息实现了图像的增强。文献[8]提出了一种改进的伽马校正函数，该方法通过图像的局部信息自适应地确定伽马函数中的参数，有效地提高了图像的对比度和亮度。

图像的梯度反映了一副图像的纹理信息，对图像的梯度场进行处理可以有效地提升图像的质量。目前针对于图像梯度场进行操作的图像增强方法其步骤主要可分为以下三步：首先获取原图像的梯度场；其次对原梯度场进行某种变换以获得目标梯度场；最后根据目标梯度场进行图像重建获得增强后的图像[9]。在文献[9]中，董丽丽、丁畅等通过对梯度场进行局部均衡化，抑制了图像高亮区域的扩散，取得了不错的增强效果。文献[10]中，赵文达等提出一种基于图像梯度直方图高斯规定化的红外图像增强算法，该算法通过对原图像的梯度直方图进行高斯函数规定化，使红外图像的成像效果更加清晰。

本文在传统伽马校正的基础上提出一种改进的伽马校正函数，该函数的参数融合了每个分块图像的全局信息。同时将改进的伽马变换作用于图像的梯度域，获取目标梯度场后再进行图像的重建。在图像重建阶段，本文通过对图像进行分块处理有效地减少了计算量，大大缩减了图像的重建时间。另外，由于本文采用分块重建策略重建了整个图像，在算法具体实施时更具实用性。

2 算法原理

2.1 梯度场的阈值分割

本文在对图像的梯度场进行增强之前首先需要对其进行梯度划分，在图像梯度场中，梯度值较小的像素点在所有的像素点中所占比例很高，体现在梯度直方图中就是左高右低的类似于“滑坡”的样式，如图1所示。针对于这种分布特点，文献[11]提出一种新的梯度阈值划分方法，即根据梯度直方图的累计分布比例来划分图像梯度场的大梯度区间和小梯度区间。本文采用这种梯度划分方法，在将图像梯度场进行划分以后，对图像的大梯度区间进行适度增强操作；对于小梯度区间的梯度进行大幅度的增强。以此来达到图像梯度场共同增强的目的。

图1 梯度直方图

2.2 改进的伽马校正算法

伽马校正早期被应用于显示器和扫描仪这一类的成像设备上，应用它的目的主要是对输入图像进行校正。伽马校正的具体形式见式(1)。

S(x,y)=C·R(x,y)γ

(1)

式中：R(x,y)表示输入图像的灰度值；S(x,y)表示经过伽马校正后的输出图像灰度值；C表示正值常数，为了将输入灰度值和输出灰度值归一化为[0,1]的范围，C值通常取1，因此式(1)可简化为S(x,y)=R(x,y)γ；式中的γ被称为伽马校正系数，随着其取值的不同图像会产生不同的校正效果。

从图2可以看出，当γ>1时，函数输出值会随输入值的增大而减小，具体表现为图像的灰度值被压缩，图像变得更暗；反之，当γ<1时，伽马函数会增大原图像的灰度值，图像会变得更亮；当γ=1时，伽马函数将退化为恒等变换，不会对输入图像产生作用。因此为了达到增强梯度场的目的，本文所提出的改进的伽马函数主要在γ<1的系数范围内对原图像梯度场进行操作。

图2 传统伽马校正函数图像

传统伽马变换对于图像的梯度场具有一定的增强作用，但其在加深梯度图像背景的同时可能会拉低图像中目标物的梯度值[12]。此外，在传统伽马变换中，通常使用单个固定的伽马校正系数对整张图像进行处理，这样只能对目标图像进行局部的增强。因此针对于以上两点局限性，本文提出一种基于图像梯度场局部信息的伽马变换函数，见式(2)。

(2)

其中：

(3)

式中：Gout(x,y)表示伽马函数输出的图像梯度值；Gin(x,y)表示输入的图像梯度值，在本文算法中，输入梯度必须要归一化至[0,1]的区间范围；Gmax表示原图像的梯度最大值；T为2.1节所述的梯度图像分割的梯度阈值，其数值表示的是梯度值而非梯度图像分割的比例；θ和μ分别代表原图像的梯度标准差和梯度平均值；β1和β2是调整系数。

在式(2)和式(3)中，β1和β2需要进行设定，其余参数来源于梯度图像的统计信息。图3是基于图像统计数据绘制的函数图像。

图3 改进的伽马函数图像

从图3可知，在不同的梯度分割区间内，本文提出的伽马函数对原图像的梯度值采用了不同的放大策略。与大梯度值变换曲线相比，小梯度值变换曲线要更加陡峭一些。这体现出改进的伽马函数对大梯度区间进行了一定程度的放大，而对于小梯度区间产生了更强的放大效果。这与2.1节中所提出的梯度分区间变换思路相符合。需要注意的是，由于大梯度区间的梯度值比较大，在经过函数处理后会出现输出值过大的情况，因此本文在算法实现阶段对其输出值进行了限制，见式(4)。

Gout(x,y)=η·Gmax,Gout(x,y)>Gmax

(4)

对于低照度图像，η取值一般在1～2.5范围内，本文取2。

2.3 图像重建

对原图像梯度场进行增强处理后会得到目标的梯度场G，因此需要寻找一个梯度场与目标梯度场 最佳近似的灰度图像U[13]。在数学上该过程等价于：在所有二维函数空间中，寻找一个梯度在最小二乘意义下与目标梯度场G(x,y)最接近的函数U(x,y)，即求式(5)的泛函最小值[14]：

(5)

利用变分法[15]中的Euler-Lagrange方程对式(5)进行化简整理最终可得式(6)：

(6)

可以看出式(6)就是泊松方程，因此对目标梯度场的重建问题可以转化为求泊松方程的数值解问题。而对泊松方程进行数值求解实际上就是求解以下的线性方程组：

LU=divG

(7)

式中：L是Laplacian系数矩阵，该矩阵可以通过卷积方式构建出来[16]；U是所求图像矩阵的列向量形式。G是目标梯度场的散度向量，求解式见式(8)。

divG=Gx(x,y)-G(x-1,y)+Gy(x,y)

-G(x,y-1)

(8)

上述方法虽然在算法的实现方面较为简单，但是仍然存在较大的问题。以本文采用的512×512的图像为例，在进行图像复原时，需要将目标梯度场的散度矩阵按行拉伸成(262144,1)的列向量，因此与之相匹配的系数矩阵L的阶数会迅速增加到(262144,262144)。过大的系数矩阵会造成数据存储的困难，而且也很难对其进行稀疏化。另外，在线性方程组的计算方面也会消耗大量的计算资源[10]。针对以上问题，文献[5]提出一种泊松方程快速解法，该方法是将二维的Laplacian系数矩阵分解为两个正交方向上的1维矩阵，并对其进行正弦变换，最后再利用Kronecker直积方法，将两个1维向量组合为2维Laplacian系数矩阵。该方法解决了2维Laplacian系数矩阵的生成问题，但系数矩阵的阶数并未下降。在文献[11]中，丁畅、董丽丽等在文献[5]的基础上提出一种矩阵变换法，该方法通过三角矩阵的相似变换对角化，大幅度地降低了系数矩阵的阶数，计算速度显著提高，但该方法在工程实现时需要多次矩阵变换稍显复杂，因此希望能采用较为简单的策略解决这个问题。

本文采用对图像进行分块重建的方式降低计算量。以512×512的图像为例，如果将原图像划分为4×4的图像块，则单块图像的系数矩阵会迅速下降为(16,16)，不但解决了系数矩阵存储的问题，还大大地缩减了图像的重建时间。另外，这样做避免了算法实现过程中的复杂变换，更具实际过程意义。

2.4 算法的总体流程

算法的总体流程见图4。本文采取的方法在具体实现时存在以下几点说明：

图4 算法的总体流程

(1)原图像的梯度需要分成水平和竖直两个方向分别求出。另外，因为图像的梯度值可能存在正负值及零值，因此需要用两个矩阵分别保存两幅梯度图像的正负信息。

(2)在求取散度之前，需要将原梯度图像的方向信息进行还原。

(3)求取图像的散度信息以后，需要将其按行压缩为列向量，再通过2.3节方法重建后需要将图像按照分块的尺寸重构成分块图像。

3 实验评价

3.1 本文方法效果的主观评价

在图像处理领域，对图像具有增强作用的常用算法有直方图均衡化和基于“视网膜大脑皮层理论”的Retinex算法(SSR、MSR、MSRCR)。因此本文在实验中选用文献[5]算法和以上两种算法与本文方法进行比较。通过比较图5至图7发现，本文方法可以比较有效地提升MRI的对比度和细节轮廓。另外，直方图均衡化等方法在处理图像时对背景噪声进行了增强。相比较而言，本文方法噪声幅值更小，在视觉体验上更加“干净”。

3.2 本文方法效果的客观评价

目前常用的图像质量评价方法有很多，本文采用信息熵[17]、图像标准差和峰值信噪比作为算法的评价指标。信息熵一般用来评价图像中所含信息的丰富程度，表1是图5至图7中各类算法的信息熵比较结果。可以看出，本文方法在一定程度上提升了原图像的信息熵。由于直方图均衡化等算法在图像增强的过程中增大了图像背景噪声的强度，因此其较高的信息熵实际上体现了其过度增强的本质。

表1 不同算法的信息熵对比

图像的标准差是描述图像细节的指标，标准差值越大说明图像的细节表现得越清楚。表2的数据结果显示了本文方法对于MRI具有较好的增强效果。

表2 不同算法的标准差对比

为了衡量在图像增强过程中各类算法是否同步地对噪声进行了增强，本文采用峰值信噪比(PSNR)来进行各类算法结果的评价。峰值信噪比在图像处理领域中是一种非常常用的客观评价指标，它常用于评价图像处理结果的失真程度。峰值信噪比越大，图像处理结果失真越小，与原图像的相似性越大。从表3中的数据可以看出，本文方法的处理结果与原图像相似性更大，说明本文方法只是对图像特征进行增强，而并未进行图像噪声增强。

表3 不同算法的PSNR对比

3.3 本文算法速度提升评价

由2.3节可知，本文通过小分块重建的方式对增强后的梯度场进行图像重建。在实验中发现，采用分块的策略可以较大幅度地缩短算法的执行时间，同时也能降低算法对内存等计算资源的消耗。在图5(e)、图6(e)和图7(e)中，本文采用4×4的分块尺寸进行重建，程序耗时分别是：5.24 s、5.38 s和5.49 s。由于本文采用的512×512分辨率，图像尺寸过大，用传统方法很难对其进行整体重建，因此，另外选择了一幅200×200分辨率的图像来进行对比测试。具体的计算时长参照表4。

图5 膝盖部位的MRI增强效果图及直方图

图6 脑部MRI增强效果图及直方图

从表4中可以看出，伴随着分块尺寸的缩小，本文算法的执行速度得到较大的提升。这说明本文采用的分块重建的策略是成功的。在算法速度上基本满足了实时性的要求。

表4 不同分块尺寸的时长对比(单位：s)

3.4 方法的不足之处

(1)本文方法在对梯度场进行改进的伽马函数处理时未能完全实现参数自适应，因此针对不同种类的图像需要调整不同的β1,β2参数值，但是过大的β1,β2可能会导致大、小梯度区间的变换曲线置于恒等变换曲线以下，即会对图像梯度产生压缩效果，因此针对于MRI图像，如果参数取以下的范围值可以获得较好的效果：β1=6～9,β2=β1/2.5。

(2)通过观察图5至图7发现，在采用本文方法进行图像分块合并以后，图像效果基本符合人眼的视觉要求。但如果在处理过程中分块尺寸过大，经过处理以后的图像块之间的灰度值会产生较大的差异，合并以后会产生比较明显的块状效应，这个弊端是本文采用的分块策略导致的。

以上不足之处有待后续改进。

4 结论

针对部分MRI图像中存在的细节模糊不清、对比度低的问题，本文提出一种基于改进伽马校正的图像增强算法。该算法首先在原图像梯度场上进行改进的伽马校正函数处理，然后采用小分块重建的方式对增强后的梯度场进行重建，最后再合并小的图像块获取增强后的图像。通过与其他算法进行比较可以看出，本文提出的方法对于MRI图像具有较好的增强效果，且在运算速度上也超过了传统算法，对于本文采用的512×512分辨率图片，在i3-10110U CPU和12 G内存的硬件条件下，本文算法最快的处理时间约为5.24 s。此外，本算法也为传统梯度图像重建算法的改进提供了一种新的思路。针对于本文方法的不足之处，有待后续进行持续的改进，以使其在具体应用方面更加成熟。