MRI图像增强和灰度插值

题    目： MRI图像增强和灰度插值 

专    业：     生物医学工程       

班    级：      生医2班           

学    号：              

姓    名：                  

指导教师：                 

天津医科大学 生物医学工程学院

年     月    日

一. 开题背景

  在社会生活和科研生产中, 人们随时随地都要接触图像。图像信息是人类认识世界的重要信息来源。而图像处理就是一门关于图像的综合学科, 它汇聚了光学、 电子学、 数学、 摄影技术和计算机技术等众多学科【1】。

    所谓图像处理,就是通过某些数算对图像信息进行加工, 以满足人的视觉心理和实际生活需要。医学图像处理作为一门具有很强实用价值和广泛应用前景的交叉学科,吸引了电子、数学、物理等诸多不同专业人员从事于这一热点领域的研究。医学图像具有信息量大、图像模糊、处理困难等特点,借助于图形图像技术的快速发展和应用,医学图像的成像质量和显示方法得到了极大的改善,大大提高临床诊断的准确性和正确性,这样既能充分挥了医学影像设备的效力和潜能,又能促进了诊疗水平的提高,改善了人们群众的健康况。因此对医学图像处理方法的研究具有十分重要的理论意义和现实意义【2】。

    磁共振成像MRI可以很好地识别大脑的灰质和白质密度接近的软组织，从而深受医生的青睐。其临床应用越来越广，然而受其他因素影响。原始的MRI图像可能出现值脉冲噪声、伪影、图像模糊等。因此需要对MRI图像进行增强处理，凸显重要的医学特征改善图像的质量，以方便医生对病情做出准确的判断。

二．课题目的

在学习了数字信号处理和医学图像处理等课程的基础上，以MRI图像为例，利用Matlab软件，学习医学图像的读取保存等基本操作，图像的灰度映射增强(直方图均衡化、直方图规定化和Gamma校正)以及常用的三种灰度插值算法：最近邻算法，双线性插值法，三次多项式插值。

三．课题研究的主要内容

1、原始图像是MRI脑图像，分别为t1加权像，t2加权像和pd质子加权像。了解所给MRI图像的图像信息，如图像宽度、高度和文件类型等；

2、根据图像的不同类型，把原始MRI图像转为灰度图像，将图像拆分为大小相等的t1, t2, pd三幅图像并保存；

3、绘出t1, t2, pd三幅图像的灰度直方图，分析图像的灰度分布情况，并分别利用直方图均衡化、直方图规定化和Gamma校正对三幅图像进行增强处理，并对结果加以比较； 

4、对t1, t2, pd三幅图像每三点取一点进行采样，然后分别利用最近邻法，双线性法及三次多项式插值恢复图象，并对这三种插值方法结果进行比较。

四．原理和方法

1、利用Matlab软件的imread函数读入图像，再imfinfo函数获取所给MRI图像信息；

2、根据得到的图像的矩阵信息将图像拆分为大小相等的t1, t2, pd三幅图像并保存；

3、利用imhist函数绘出三幅图像的灰度直方图，再分别用histeq函数，自定义的                       

   twomodegauss函数，imadjust函数对三幅图像进行直方图均衡化，直方图规定化，                           

Gamma 校正增强处理；

4、对t1, t2, pd三幅图像每三点取一点进行采样，再利用imresize函数进行最近邻法，双   

线性法及三次多项式插值恢复图象【3】。

五．步骤与结果

1、利用imfinfo函数获取所给MRI图像信息，程序见附录（1），结果如下:

 Width: 520

 Height: 284

 ColorType: 'indexed'

由此可知原始MRI图像的宽为520列，高为284行，图像类型是索引图像；

2、利用ind2gray函数将原始图像转换为灰度图像，程序见附录（1），结果如图1:

 图1  原始MRI的灰度图像

3、根据得到的图像的矩阵信息将图像拆分为大小相等的t1, t2, pd三幅图像，程序见附

 录（1），结果如图2.1~2.3:

 图2.1  原始t1加权像的灰度图像

 图2.2  原始t2加权像的灰度图像

 图2.3  原始pd质子加权像的灰度图像                    

4、绘出t1, t2, pd三幅图像的灰度直方图，程序见附录（2），结果如图3.1~3.3：

 图3.1  t1加权像的灰度直方图  

 图3.2  t2加权像的灰度直方图

 图3.3  pd质子加权像的灰度直方图

结果分析：由三幅加权像的灰度直方图可以看出，脑组织的灰度值主要分布在40-150  

 之间。

5、图像增强处理

（1）直方图均衡化

     利用histeq函数对三幅加权像的灰度图像进行均衡化处理，程序见附录（3），结果如图4.1~4.6：

 图4.1  t1加权像均衡化后图像

 图4.2  t1加权像均衡化后图像的灰度直方图

 图4.3  t2加权像均衡化后图像

 图4.4  t2加权像均衡化后图像的灰度直方图

 图4.5  pd加权像均衡化后图像

 图4.6  pd加权像均衡化后图像的灰度直方图

结果分析：

    直方图均衡化的基本思想是把原始图像的直方图变换成均匀分布的形式，这样就增加了图像灰度值的动态范围，从而达到了增强图像整体对比度的效果。但这种方法的缺点是它对处理的数据不加选择，它可能会增加背景杂讯的对比度，并且降低有用信号的对比度，变换后图像的灰度级减少，某些细节消失。从而导致三幅加权像均衡化效果不佳。

（2）直方图规定化

     利用自定义的一个双峰高斯函数作为规定化函数对三幅加权像的灰度图像进行规定化处理，双峰高斯函数如图5.1，程序见附录（4），规定化程序见附录（5），结果如图5.2~5.7：

 图5.1  双峰高斯函数

 图5.2  t1加权像规定化后的图像

 图5加权像规定化后的图像的灰度直方图

 图5.4  t2加权像规定化后的图像

 图5加权像规定化后的图像的灰度直方图

 图5.6  pd加权像规定化后的图像

 图5加权像规定化后的图像的灰度直方图

结果分析：

    直方图规定化就是通过一个灰度映射函数，将原灰度直方图改造成所希望的直方图。所以，直方图修正的关键就是灰度映像函数。这里选择双峰高斯函数作为灰度映射函数，规定化后图像层次清晰，效果较好。

(3) Gamma 校正

   利用imadjust函数对三幅加权像的灰度图像进行Gamma 校正，这里以t1加权像的结果为例，程序见附录（5），结果如图6.1~6.8：

 图6.1  gamma=0.2时t1加权像的校正图像               

 图6.2  gamma=0.4时t1加权像的校正图像

 图6.3  gamma=0.6时t1加权像的校正图像

 图6.4  gamma=0.8时t1加权像的校正图像

 图6.5  gamma=1.4时t1加权像的校正图像

 图6.6  gamma=1.8时t1加权像的校正图像

 图6.7  gamma=2.2时t1加权像的校正图像

 图6.8  gamma=2.4时t1加权像的校正图像

结果分析：

    当gamma>1时，输入中较宽的低灰度范围被映射到输出中较窄的灰度范围，输入中较窄的高灰度范围被映射到输出中较宽的灰度范围；这里gamma=1.4时校正图像效果较好。

6、图像插值

（1）对t1, t2, pd三幅加权图像每三点取一点进行采样，然后利用imresize函数分别进行最近邻法，双线性法及三次多项式插值恢复图象，比较不同放大倍数时不同插值方法的效果，这里以t1加权像的插值结果为例，程序见附录（6），结果如图7.1~7.9：

 图7.1  放大倍数为1时最近邻插值图像

 图7.2  放大倍数为1时双线性法插值图像

 图7.3  放大倍数为1时三次多项式插值图像

 图7.4  放大倍数为2时最近邻插值图像

 图7.5  放大倍数为2时双线性插值图像

 图7.6  放大倍数为2时三次多项式插值图像

 图7.7  放大倍数为3时最近邻插值图像

 图7.8  放大倍数为3时双线性插值图像

 图7.9  放大倍数为3时三次多项式插值图像

结果分析：

         放大倍数较小时，三种插值方法效果差别不大；当图像放大倍数较大时，双    

         线形插值法和三次多项式插值效果较好，最近邻插值效果较差。

(2)插值图像和原始图象进行相减，比较两幅图象在灰度上的差别，统计差值图像的平均  值和方差等指标，对插值图象的好坏进行评估：

  插值图像与原始图像相减，这里以t1加权像为例，程序见附录（7），结果如图8.1~8.6：

 图8.1  t1加权像、最近邻插值图像、相减后图像

 图8.2  最近邻插值图像与t1加权像相减后图像的灰度直方图

 图8.3  t1加权像、双线性插值图像、相减后图像

 图8.4  双线性插值图像与t1加权像相减后图像的灰度直方图

 图8.5  t1加权像、三次多项式插值图像、相减后图像

 图8.6  三次多项式插值图像与t1加权像相减后图像的灰度直方图

结果分析：

    由相减图像的灰度直方图可以看出，最近邻插值法可能会造成插值生成的图像在灰度上的不连续，在灰度变化的地方可能出现明显的锯齿状。

    双线性插值法与三次多项式插值法生成的图像在灰度没有不连续的缺点，灰度变化的地方较平滑，结果基本令人满意。

差值图像的平均值和方差，程序见附录（8)，结果如表1~表3：

               表1  插值图像与t1加权像差值图像的均值和方差

               表3  插值图像与pd加权像差值图像的均值和方差

    由表1~表3可以看出差值图像的均值与方差都是三次多项式插值法最小，双线性插值法次之，最近邻插值法最大，说明三次多项式插值图像效果最好，最近邻插值图像效果最差。

插值图像与原始图像的相关系数，程序见附录（9），结果如表4~表6：

                   表4  插值图像与t1加权像的相关系数

    由表4~表6可以看出插值图像与原始图像的相关系数都大于0.9，且相差很小，此指标不适合评估不同插值方法的插值效果。

三种不同插值算法的运行时间，程序见附录（10），结果如表7：

               表7  三幅加权像三种插值法循环100次运行时间

    由表7可以看出三次多项式插值计算成本最高，双线性插值次之，最近邻插值计算成本最小。

六．结论

1.直方图均衡化可能会增加背景杂讯的对比度，并且降低有用信号的对比度，变换后图 像的灰度级减少，某些细节消失。

2.放大倍数较小时，最近邻、双线性、三次多项式三种插值方法效果差别不大；当图像放大倍数较大时，双线性插值法和三次多项式插值效果较好，最近邻插值效果较差。

3.最近邻插值法会造成插值生成的图像在灰度上的不连续，在灰度变化的地方可能出现明显的锯齿状。双线性插值法与三次多项式插值法生成的图像在灰度没有不连续的缺点，灰度变化的地方较平滑，结果基本令人满意。

4.从差值图像的均值与方差两个指标来看，都是三次多项式插值法最小，双线性插值法次之，最近邻插值法最大，说明三次多项式插值图像效果最好，最近邻插值图像效果最差。

5.三次多项式插值计算成本最高，双线性插值次之，最近邻插值计算成本最小。

七．讨论

图像是离散函数，同时近似运算存在误差，规定化变换只能接近参考直方图，不可能完全相同。现有的图像增强的方法各有利弊，同一种方法能突出图像的某些特征，又会掩盖或消除图像别的特征，因此应根据需要加以选用。

参考文献

【1】罗笑南, 王若梅. 计算机图像学. 广州: 中山大学出版社, 1996

【2】239000 石永华 . 基于contourlet变换的 MRI医学图像增强 . 滁州学院学报，2011年10月

【3】章毓晋.图像处理.第三版.清华大学出版社.北京

附录

（1）

s0=size(f);

f1=ind2gray(f,map);figure;imshow(f1);title('gray image');

f11=f1(1:s0(1),1:s0(2)/3);figure;imshow(f11);title('t1');

f12=f1(1:s0(1),s0(2)/3+1:s0(2)*2/3);figure;imshow(f12);title('t2');

f13=f1(1:s0(1),s0(2)*2/3+1:s0(2));figure;imshow(f13);title('pd');

（2）

figure(1);imhist(f11);

figure(2);imhist(f12);

figure(3);imhist(f13);

（3）

g1=histeq(f11);figure(1);imshow(g1);title('t1');

figure(2);imhist(g1);

g2=histeq(f12);figure(3);imshow(g2);title('t2');

figure(4);imhist(g2);

g3=histeq(f13);figure(5);imshow(g3);title('pd');

figure(6);imhist(g3);

（4）

function p=twomodegauss(m1,sig1,m2,sig2,A1,A2,k)

c1=A1*(1/((2*pi)^0.5)*sig1);

k1=2*(sig1^2);

c2=A2*(1/((2*pi)^0.5)*sig2);

k2=2*(sig2^2);

z=linspace(0,1,256);

p=k+c1*exp(-((z-m1).^2)./k1)+c2*exp(-((z-m2).^2)./k2);

p=p./sum(p(:));

p1=twomodegauss(0.15,0.05,0.75,0.05,1,0.07,0.002);

figure;

plot(p1);

xlabel('灰度值');

ylabel('概率');

（5）

g11=imadjust(f11,[],[],0.2);

g12=imadjust(f11,[],[],0.4);

g13=imadjust(f11,[],[],0.6);

g14=imadjust(f11,[],[],0.8);

g15=imadjust(f11,[],[],1.4);

g16=imadjust(f11,[],[],1.8);

g17=imadjust(f11,[],[],2.2);

g18=imadjust(f11,[],[],2.4);

figure(1);imshow(g11);title('gamma=0.2');

figure(2);imshow(g12);title('gamma=0.4');

figure(3);imshow(g13);title('gamma=0.6');

figure(4);imshow(g14);title('gamma=0.8');

figure(5);imshow(g15);title('gamma=1.4');

figure(6);imshow(g16);title('gamma=1.8');

figure(7);imshow(g17);title('gamma=2.2');

figure(8);imshow(g18);title('gamma=2.4');

（6）

s=size(f11);m=s(1);n=s(2);

r1=f11(1:2:m,1:2:n);

r11=imresize(r1,1,'nearest');

figure(1);imshow(r11);title('nearest');

r12=imresize(r1,1,'bilinear');

figure(2);imshow(r12);title('bilinear');

r13=imresize(r1,1,'bicubic');

figure(3);imshow(r13);title('bicubic');

(7)

s=size(f11);m=s(1);n=s(2);

r1=f11(1:2:m,1:2:n);

r11=imresize(r1,[m n],'nearest');

h11=r11-f11;

figure(1);subplot(131);imshow(f11);title('t1 original');

subplot(132);imshow(r11);title('nearest');

subplot(133);imshow(h11);title('subtraction');

r12=imresize(r1,[m n],'bilinear');

h12=r12-f11;

figure(2);subplot(131);imshow(f11);title('t1 original');

subplot(132);imshow(r12);title('bilinear');

subplot(133);imshow(h12);title('subtraction');

r13=imresize(r1,[m n],'bicubic');

h13=r13-f11;

figure(3);subplot(131);imshow(f11);title('t1 original');

subplot(132);imshow(r13);title('bicubic');

subplot(133);imshow(h13);title('subtraction');

(8)

s=size(f11);m=s(1);n=s(2);

r1=f11(1:2:m,1:2:n);r2=f12(1:2:m,1:2:n);r3=f13(1:2:m,1:2:n);

r11=imresize(r1,[m n],'nearest');

r12=imresize(r1,[m n],'bilinear');

r13=imresize(r1,[m n],'bicubic');

h11=r11-f11;

average11=mean(mean(h11))

var(h11(:))

h12=r12-f11;

average12=mean(mean(h12))

var(h12(:))

h13=r13-f11;

average13=mean(mean(h13))

var(h13(:))

r21=imresize(r2,[m n],'nearest');

r22=imresize(r2,[m n],'bilinear');

r23=imresize(r2,[m n],'bicubic');

h21=r21-f12;

average21=mean(mean(h21))

var(h21(:))

h22=r22-f12;

average22=mean(mean(h22))

var(h22(:))

h23=r23-f12;

average23=mean(mean(h23))

var(h23(:))

r31=imresize(r3,[m n],'nearest');

r32=imresize(r3,[m n],'bilinear');

r33=imresize(r3,[m n],'bicubic');

h31=r31-f13;

average31=mean(mean(h31))

var(h31(:))

h32=r32-f13;

average32=mean(mean(h32))

var(h32(:))

h33=r33-f13;

average33=mean(mean(h33))

var(h33(:))

（9）

corr2(r11,f11)

corr2(r12,f11)

corr2(r13,f11)

corr2(r21,f12)

corr2(r22,f12)

corr2(r23,f12)

corr2(r31,f13)

corr2(r32,f13)

corr2(r33,f13)

（10）

tic;

for i=1:100

r11=imresize(r1,[m n],'nearest');

r21=imresize(r2,[m n],'nearest');

r31=imresize(r3,[m n],'nearest');

end

toc;

tic;

for i=1:100

r12=imresize(r1,[m n],'bilinear');

r22=imresize(r2,[m n],'bilinear');

r32=imresize(r3,[m n],'bilinear');

end

toc;

tic;

for i=1:100

r13=imresize(r1,[m n],'bicubic');

r23=imresize(r2,[m n],'bicubic');

r33=imresize(r3,[m n],'bicubic');

end

toc;