基于MATLAB的MRI及CT图像增强处理

1、共17页第 页共17页第 页课程设计说明书题目：基于MATLAB的MRI及CT图像增强处理学院（系）：年级专业：学号:学生姓名：指导教师：教师职称：目录一、课程设计目的4TOC o 1-5 h z HYPERLINK l bookmark4 二、课程设计要求4 HYPERLINK l bookmark6 三、课程设计内容4四、题目分析4 HYPERLINK l bookmark10 五、总体设计巧 HYPERLINK l bookmark12 六、程序设计51、灰度转换-52、灰度变换63、空余滤波74、频域滤波85、伪彩色增强11 HYPERLINK l bookmark14 七、心得体会1

2、2八、附录13附：燕山大学课程设计评审意见表共17页第 页共17页第 页一、课程设计目的学习使用matlab软件，并使用matlab软件编写相应的程序完成对图像的增强处理。二、课程设计的要求1熟悉和掌握MATLAB程序设计方法2掌握matlab软件在图像增强技术的处理方法3编写相应的matlab程序进行图像处理三、课程设计的内容学习MATLAB程序设计，利用MATLAB图像处理工具箱，设计和实现自己的matlab图像增强设计要求。要求：按照软件工程方法，根据需求进行程序编程，给出设计详细说明。然后按照自己拟定的功能要求进行程序设计和调试。以下几点是程序必须实现的功能。1）图像的读取和保存。2）

3、能够采用各种方法对图像进行处理，显示和对比变换前后的图像。3）编写程序并进行相应的程序说明。比较各种方法的好坏。4）图像直方图统计和直方图均衡，要求显示直方图统计，比较直方图均衡后的效果。5）能对图像加入各种噪声，并通过几种滤波算法实现去噪并显示结果。比较去噪效果。四、题目分析信息化社会中，计算机在各种信息处理中发挥着重要的作用。我们可以借助计算机，对数字图像进行处理，以达到不同的效果。根据题目的要求，除了实现要求的功能外，还有很多的功能需要用到。（1）、由于医学图像采集到的都是灰色图像，而有许多图像格式并不符合要求，因此要通过编程将一个RGB图像转换为灰度图像。（2）图像增强可以从以下几方面

4、入手：灰色增强、空运增强、频域增强、彩色增强。（3）必须对图像读取和显示这种基本处理以及相应的对比。（4）、选择MATLAB计算编程可视化，功能强大，学、易用；工具箱中包括图像处理函数,涵盖了几乎所有的图像处理技术方法;具有开放性,用户可以根据需要书写自己的函数，以满足特定的需要;能够准确读取DICOM格式文件,并有丰富的处理DICOM格式图片的各种函数。医学影像设备产生的文件为DICOM格式,普通图像处理软件无法读取此格式，若将其转换为JPG或BMP又会造成某些信息的丢失,MATLAB可以克服这一缺点。、分析一个图像的频谱特征，利用傅里叶变换，将图像从空间域变换到频域，然后进行各种处理，经过

5、高通滤波器或是低通滤波器。、由于医学图像采集到的都是灰色图像，可以通过matlab进行伪彩色增强得到相应的彩色图像。、五、总体设计本课程设计主要从以下四个方面对图像进行处理：灰度增强、空域增强、频域增强、彩色增强、1、灰度增强：灰度变换的方法包括直方图灰度变换、直方图均衡化、直方图规定化、本课程设计主要采用直方图均衡化完成图像增强处理、2、空域增强：本课程设计采用平滑滤波器来完成空域滤波，其中采用均平滤波器filter2、中值滤波器medfilt2、维纳滤波器wiener、3、频域滤波：在医学图像里采用低通滤波器进行滤波，采用高通滤波器进行边缘检出、4、彩色增强：由于MRI图像得到的都是灰色图

6、像，而人眼对彩色敏感度更高，因此可以采用伪彩色增强的方法对图像进行增强处理以得到彩色图像六、程序设计1、灰度转换由于RGB图像是三维图像，所以图像数据是一个三维数组，为了显示灰度图像，把三维图像降为二维，可以只取其中的二维数据、我在设计中采用rgb2gray函数来完成数据的转换，但是其必须指出其开始是否灰色图像，否则容易出错、具体程序如下：I=imread(naogan.jpg);ifisrgb(I)I2=rgb2gray(I);elsemsgbox(该图像是灰色图像，无需转换)；endfigure(1),imshow(I),title(原图)figure(2),imshow(I2),titl

7、e(处理后)处理后的图像如下图一所示：2、灰度变换利用图像灰度级的分布可以看出图像灰度分布的特性。如果大部分像素集中在低灰度区则图像呈现暗特性,反之则呈现亮特性。灰度变换的目的是通过改善直方图的灰度分布特性,进而改善图像的质量。灰度变换的方法包括直方图灰度变换、直方图均衡化、直方图规定化。本研究以直方图均衡化为例说明该模块的设计功能。直方图均衡化的基本思想是把原始图像的直方图变换成均匀分布的形式，这样就增加了图像灰度值的动态范围，从而达到了增强图像整体对比度的效果。均匀量化的自然图像的灰度直方图通常在低值灰度区间上频率较大，使得图像中较暗区域中的细节常常看不清楚。为了使图像清晰，可将图像的灰度

8、范围拉开，让灰度频率较小的灰度级变大，即让灰度直方图在较大的动态范国内趋于一致。理想情况下：使一输入图象转换为在每一灰度级上都有相同的像素点数(即输出的直方图是平的)。具体方法是：列出原始图像的灰度级Sk,k=O,l,，L1；其中L是灰度级的个数。统计原始图像各灰度级的像素数目nk。计算原始图像直方图各灰度级的频率数。计算原始图像的累计直方图。(5)取整计算,tk=int(N1)tkkN确定映像关系，Sknk统计新直方图各个灰度级的像素数目nk。计算新的直方图，Pt(tk)=nk/No如要处理的mri图像，通过直方图分析,原图像大面积为暗色，并且层次不清,经过直方图均衡化后,直方图的灰度间隔被

9、拉大,显得较为“平坦”,灰度层次等级增加,然后用此均衡得到理想的图像。其程序代码如下：I=imread(mri.tif);J=histeq(I);subplot(121);imshow(I);title(原图)subplot(122);imshow(J);title(直方图均衡化后);figure,subplot(121);imhist(I,64);title(原始图像直方图);subplot(122);imhist(J,64)title(均衡变换后直方图)其图像处理结果如下：*Figure1均衡化处理前后脑部mri图片均衡化前后的脑膜炎图片直方图比较3、空域增强空域滤波在处理效果上来分，可以

10、分为平滑滤波器和锐化滤波器。下面以平滑滤波器为例介绍课件功能及实现。平滑滤波器的目的在于消除混在图像中的干扰常使用的滤波器有均平滤波器filter2()、中值滤波器medfilt2()、维纳滤波器wiener()。由于医学成像设备CT机自身的原因，在CT图像中容易产生高斯白噪声，严重影像了图像的质量，为了去除此类干扰可采用维纳滤波器;如果影像本身带有盐椒噪声，可采用中值滤波器。通过给图像增加盐椒噪声(SALTPEPPET)或高斯噪(GUASSIAN)来模拟医院的ct及存在噪声的mri图像，然后采用不同滤波器进行处理对比其效果。其对应的的matlab程序如下：共17页第 页共17页第 页I=im

11、read(mri.tif);J=imnoise(I，salt&pepper,0.02);J=imnoise(I,gaussian,0,0.005);K=filter2(fspecial(average,3),I)/255;L=medfilt2(I,33);M=wiener2(I,55);%添加盐椒噪声%添加高斯噪声%均平过滤器%中值过滤器%维纳滤波器%显示图像经过对比得知经中值滤波得到的图像效果最佳。subplot(l,5,l),imshow(I),title(原图)subplot(l,5,2),imshow(J),title(添加噪声后的图片)subplot(1,5,3),imshow(K)

12、,title(均平过滤后的图片)subplot(1,5,4),imshow(L),title(中值过滤的图片)subplot(1,5,5),imshow(M),title(维纳滤波后的图片)4、频域滤波由于图像的边缘、跳跃部分及噪声颗粒是图像信号的高频分量大面积的背景区域代表了图像信号的低频分量。因此，在频域内采用低通滤波器可以进行滤波，使用高通滤波器可以进行边缘检出。其中，常用低通滤波器有理想低通滤波器、Butterworth低通滤波器、指数低通滤波器和梯形低通滤波器等;常用高通滤波器有理想高通滤波器、Butterworth高通滤波器、指数高通滤波器和梯形低通滤波器等。本次课程设计滤波器我采

13、用二阶巴特沃斯滤波器进行滤波。为了得到图像的频谱图，先要对数据进行傅里叶变换，用fft2函数对二维数据进行快速傅里叶变换，同时为了更好的观察频谱图，需要把fft2变换后的数据进行平移，利用fftshift函数，把快速傅里叶变换的DC组件移到光谱中心。这样图像能量的低频成分将集中到频谱中心，图像上的边缘、线条细节信息等高频成分将分散在图像频谱的边缘。（1）咼通滤波器通过高通滤波器进行边缘检出，其程序代码如下：subplot(l,2,l),imshow(I),title(原图)subplot(l,2,2),imshow(K);title(高通滤波后的图像)I=imread(mri.tif);f=d

14、ouble(I);k=fft2(f);g=fftshift(k);M,N=size(g);nn=2;d0=4;m=fix(M/2);n=fix(N/2);fori=1:Mforj=1:Nd=sqrt(i-m)A2+(j-n)A2);ifd=d0h=0;elseh=1;endresult(i,j)=h*g(i,j);endendresult=ifftshift(result);J=ifft2(result);K=uint8(real(J);%读入图像：%数据类型转换%傅里叶转换%转换数据矩阵%获取图像大小%二健巴特沃斯高通滤波器%截止频率为4%计算低通滤波器的传递函数%转换数据矩阵%傅里叶反转换

15、%转换数据类型%显示图像其相应的对比图像如下:(2)通过低通滤波器米用低通滤波器可以进行滤波，其相应的图像如下:%读入图像：%数据类型转换%傅里叶转换%转换数据矩阵%获取图像大小%二健巴特沃斯低通滤波器%截止频率是10I=imread(mri.tif);f=double(I);k=fft2(f);g=fftshift(k);M,N=size(g);nn=2;d0=10;m=fix(M/2);n=fix(N/2);fori=1:Mforj=1:N%计算低通滤波器的传递函数d=sqrt(i-m)A2+(j-n)A2);h=1/(1+0.414*(d/d0)人(2*nn);result(i,j)=h

16、*g(i,j);end%转换数据矩阵%傅里叶反转换%转换数据类型%显示图像endresult=ifftshift(result);J=ifft2(result);K=uint8(real(J);subplot(1,2,1),imshow(I),title(原图)subplot(1,2,2),imshow(K);title(低通滤波后的图像)其相应的对比图像如下图所示:5、伪彩色增强医学影像如CT、MR、B超、X光片都是灰度图像，尽管这些设备的成像质量很高，可以将灰度等级分成2000多阶，但人眼只能分别出其中16个灰度等级。若将2000个灰度等级,划分为16个灰阶则,每个灰阶所能分辨的CT值为2

17、000/16=125Hu。即相邻两组织CT值相差125时人眼才能将二者区分出来，若小于此数值，处于同一灰阶则不能分辨。而人体组织的CT值在相差几个HU单位时(35Hu)就有重要的诊断意义。然而，人眼对彩色的敏感程度远远高于对灰度的敏感程度，利用人眼的这一视觉特性，可以对医学图像进行伪彩色处理，使病灶部分能够较清晰的显现出来。灰度图像病灶部分模糊不清，但对其进行伪彩色处理后,可以较清晰地辨认出病灶的轮廓和大小。其主要实现方法为，首先读入灰度图像，然后将灰度图像分层,对图像数组进行等分层处理，最后利用PColor()进行彩色变换。其相应的程序代码如下：I=imread(mri.tif);subpl

18、ot(1,2,1);imshow(I);title(灰度图像)%将灰度图像分成8层%对图像数组进行等分层处理%进行彩色变换I=double；J=floor(I/32);J1=floor(J/4);J2=rem(floor(J/2),2);J3=rem(J,2);PColor(:,:,1)=J1;PColor(:,:,2)=J2;PColor(:,:,3)=J3;subplot(1,2,2);imshow(PColor);title(伪彩色增强后图像)其相应的对比图如下图所示:共17页第 页共17页第 页七、心得体会此次安排的matlab设计内容丰富，贴近专业知识，让我体会到了信息社会学会使用软

19、件处理工程问题的重要性，同时也是对专业知识的一次应用，了解了许多新的医学知识及图像处理的方法。通过这周的课程设计，不但使我熟悉了matlab课程设计的基本思想和基础知识，初步掌握了matlab软件关于图像处理方面的应用，而且更为深入的体会了电子信息在现代高科技信息产业领域中的重要地位。通过简单的软件编程，提高了独立思考能力和自学能力，增强了我的动手能力，当然在设计过程中我也遇到了很多问题，更是让我认识到谨慎态度的可贵。如开始时到导入一张三位图片，怎么做也无法调用imhist函数，在那里苦思冥想的一上午也没解决问题，最终才在老师的帮助下才明白matlab处理的是二维图像，最终采用图像转换函数rg

20、b2gray才得到解决；同时在设置滤波器时不明白怎么自己设置，最终查阅相关课本才完成这些问题在老师的耐心指导和帮助下都得到解决，最终顺利完成此次课程设计，得到了比较满意的结果。通过这样的练习，我发现自己学习的书本上的知识尚有欠缺，对芯片的了解只有通过自己亲手检验后才能加深。这样宝贵的经验还让我认识到合作的重要性，在整个过程中我得到了很多帮助包括老师的建议和同学的提醒，在此向热情的老师和帮助过我的同学表示感谢。附录课程设计的总程序ClearallClc%图像转换I=imread(naogan.jpg);ifisrgb(I)I2=rgb2gray(I);elsemsgbox(该图像是灰色图像，无需

21、转换);endfigure(l),imshow(I),title(原图)figure(2),imshow(I2),title(处理后)%灰度增强I=imread(mri.tif);J=histeq(I);subplot(l2l);imshow(I);title(原图)subplot(l22);imshow(J);title(直方图均衡化后);figure,subplot(l2l);imhist(I,64);title(原始图像直方图);subplot(l22);imhist(J,64)title(均衡变换后直方图)%空域增强%添加盐椒噪声%添加高斯噪声%均平过滤器%中值过滤器%维纳滤波器%显示

22、图像I=imread(mri.tif);J=imnoise(I,salt&pepper,0.02);J=imnoise(I,gaussian,0,0.005);K=filter2(fspecial(average,3),I)/255;L=medfilt2(I,33);M=wiener2(I,55);subplot(1,5,1),imshow(I),title(原图)subplot(1,5,2),imshow(J),title(添加噪声后的图片)subplot(1,5,3),imshow(K),title(均平过滤后的图片)subplot(1,5,4),imshow(L),title(中值过滤的

23、图片)subplot(1,5,5),imshow(M),title(维纳滤波后的图片)%频域增强I=imread(mri.tif);f=double(I);k=fft2(f);g=fftshift(k);M,N=size(g);nn=2;d0=4;m=fix(M/2);n=fix(N/2);fori=1:Mforj=1:Nd=sqrt(i-m)人2+(j-n)人2);ifd=d0%读入图像:%数据类型转换%傅里叶转换%转换数据矩阵%获取图像大小%二健巴特沃斯高通滤波器%截止频率为4%计算低通滤波器的传递函数%转换数据矩阵%傅里叶反转换%转换数据类型%显示图像%读入图像:%数据类型转换%傅里叶转换%转换数据矩阵%获取图像大小%二健巴特沃斯低通滤波器%截止频率是10%计算低通滤波器的传递函数h=0;elseh=1;endresult(i,j)=h*g(i,j);endendresult=ifftshift(result);J=ifft2(result);K=uint8(real(J);subplot(l