医学图像处理实验指导书(2013)

医学图像处理

实验指导书

李振伟

2013年4月

实验1 数字图像与医学图像基础

一、实验目的

1、掌握图像大小与分辨率的概念及其与图像清晰度的关系；

2、掌握图像直方图的概念及其特点；

3、了解颜色模型及其关系；

4、熟悉用Photoshop来了解图像的性质和特点。

二、实验原理

1．图像的大小与分辨率

图像分辨率指图像中存储的信息量，通常用“像素/英寸” ppi表示，在图像尺寸不变的情况下，高分辨率的图像比低分辨率图像包含的像素多，像素点较小，因而图像更清晰。

2．图像直方图

直方图的概念：对一幅图像所包含的全部像素的灰度做统计，并以横坐标表示灰度值，纵坐标表示图像中具有该灰度的像素个数，这样绘制出的曲线叫做图像的灰度直方图。

直方图的特点：

直方图反映图像灰度的统计性质，不包含空间位置信息。   

从直方图中可看出图像的总体性质。

3、颜色模型

RGB颜色模型

R代表红色，G代表绿色，B代表蓝色。RGB颜色模型是相加混色，该模型主要用于电脑监视器。

HSB颜色模型

Hue(色相)：从物体反射或透过物体传播的颜色。在 0 到 360 度的标准色轮上，按位置度量色相。在通常的使用中，色相由颜色名称标识，如红色、橙色或绿色。

Saturation(饱和度)：指颜色的强度或纯度。饱和度表示色相中灰色分量所占的比例，它使用从 0%（灰色）至 100%（完全饱和）的百分比来度量。在标准色轮上，饱和度从中心到边缘递增。

Brightness(亮度)： 颜色的相对明暗程度，通常用从 0%（黑色）至 100%（白色）的百分比来度量。

YIQ颜色模型

Y表示亮度，I表示橙－青色，Q表示绿－品红色。该模型主要用于电视监视器。

RGB模型和YIQ模型的关系：

                                       

 CMY颜色模型

C表示青色，M表示品红色，Y表示黄色，该模型主要用于彩色喷墨打印机。

三、实验内容与方法

1．改变图像的大小与分辨率

首先读入一幅图像，然后选择“图像/图像大小”菜单，在打开的图像大小对话框内输入数据设置像素大小与文档大小。

2． 查看图像的直方图

分别读入灰度图像和彩色图像，选择“窗口/直方图”菜单，查看它们的直方图。选择不同的通道，观察直方图的变化。

3．图像的亮度和对比度调整

读入一幅图像，选择“图像/调整/亮度/对比度”菜单，对图像的亮度与对比度进行调整，注意两者之间的相互关系，设置合适的亮度值与对比度值，使图像的视觉效果最佳。

4． 图像模式的转换

读入一幅图像，选择“图像/模式”菜单，分辨选中“索引颜色”、“RGB颜色”、“CMYK颜色”、“Lab颜色”，查看图像的颜色变化。

5． 图像曲线调整

打开一幅灰度图像（或打开一幅彩色图像，然后选择“图像/模式/灰度”将其转换成灰度图像），然后选择“图像/调整/曲线”，打开“曲线”控制面板，用鼠标拖动曲线进行调整，查看图像的变化。总结图像随曲线的变化规律。

6． 其他

分别选中“图像/调整”菜单下的子菜单，观察图像的变化，根据观察到的结果总结各菜单的作用及功能。

7 综合应用

对于给定图像pout.tif和forest.tif，分别用Photoshop进行以下处理：

查看图像的基本信息（包括大小、分辨率等）；

显示图像的直方图；

根据直方图分析这两幅图像的特点；

分别选用两种不同的方法对这两幅图像进行处理，使处理后的图像图像更加清晰，对比处理前后图像直方图的变化。（提示：参照实验内容3和5）。

四、思考题

    1．图像的大小与分辨率如何表示？像素大小和文档大小有什么异同？

2．什么是直方图？从直方图是否可以看出图像的整体特征？彩色图像的直方图是如何显示的？

3．调整图像的亮度与对比度有什么用途？

4．当分别选中不同的颜色模式时，图像颜色发生了什么变化？

5．图像曲线与图像之间存在什么关系？

五、实验要求

1、以上实验按顺序进行；

2、在进行每一实验内容时，除了熟悉操作之外，还应注意思考题，操作应围绕着思考题进行；

3、实验内容7 综合应用的实验过程与结果撰写在报告中；

4、实验报告内容还应包括以上5个思考题的回答；

5、实验报告完成。 

实验2 MATLAB医学图像基本操作

一、实验目的

1、掌握图像的插值方法及其特点；

2、掌握图像的基本几何变换处理；

3、熟悉用MATLAB对数字图像进行标准显示与特殊显示以及几何变换等基本操作。

二、实验原理

插值通常是利用曲线拟合的方法，通过离散的采样点建立一个连续函数，用这个重建的函数便可求出任意位置的函数值。

1.最近邻插值

概念：(U0,V0）点的灰度值等于离它最近点的灰度值,即： 

优点：简单，快捷。

缺点：误差大。当用这种方法实现大倍数处理时，在图像中可以明显地看出块状效应。

2.双线性插值

用[S]表示不超过S的最大整数，则

第一步： 

第二步：

第三步：

优点：插值效果较好。

缺点：计算量较大，并且插值后的图像变模糊，从而损失一些细节。

3.双三次插值

对于灰度变化规律较为复杂的图像，可以采用同一直线方向上更多采样点的灰度值对该点进行非线性插值。

步骤：

第一步：求a,b,c,d四点处的灰度值；

第二步：对a,b,c,d四点在垂直方向上做三次多项式内插，求得插值点的灰度值；

特点： 由16个邻近点求得，插值精度高，计算量大。

三、用到的MATLAB函数

1、imshow

功能：显示读取的图像

调用格式：

    imshow(I,n);  %I代表所显示的灰度图像的数据矩阵，n为整数，代表所要显示图像的灰度等级数，默认值为256；

   imshow(I,[low,high]);  %[low,high]为图像数据的值域。

2、colorbar 

功能：添加颜色条

调用格式：colorbar;  colorbar(‘vert’);  color(‘horiz’);

3.montage

功能：显示多帧图像

调用格式：

montage(I);

montage(X,map);

montage(RGB);

h=montage((…);

4.subimage

功能：在一个图形窗口中显示多幅图像

调用格式： 

subimage(X,map);

subimage(I);

subimage(RGB);

subimage(x,y, …);

h=subimage(…);

5.imresize

功能：对图像进行放大和缩小处理

调用格式：

B=imresize(A,m,method);  % A为原图句柄，m为放大倍数，method用于指定插值方法，可选的值为‘nearest neighbor’，‘bilinear’，‘bicubic’，默认值为‘nearest’；

B=imresize(A,[mrows,mcols],method);  ％返回一个mrows行、ncols列的图像，若与原图的长宽比不同，则图像会产生变形；

6.imrotate

功能：对图像进行旋转处理

调用格式：

B=imrotate(A,angle,method);  %method含义同上

B=imrotate(A,angle,method,’crop’);  %一般来说，旋转后的图像会比原图大，用户可以指定‘crop’参数对旋转后的图像进行剪切（取图像的中间部分），使返回的图像与原图像大小相同。

7.imcrop

功能：对图像进行剪切处理。用函数imcrop实现对图像的剪切。该操作剪切的是图像中的一个矩形子图，用户可以通过参数指定这个矩形四个顶点的坐标，也可以交互地用鼠标选取这个矩形。

调用格式：

I2=imcrop(I,rect);  

X2=imcrop(X,map,rect);

RGB2=imcrop(RGB,rect);  %按指定的矩形框剪切图像

四、实验内容与方法

（一）验证性实验

1.显示2帧图像

A1=imread(‘rice.tif’);

A2=imread(‘testpat1.tif’);

A=cat(2,A1,A2);

montage(A);

2．在一个图形窗口中同时显示4幅图像

I1=imread(‘rice.tif’);

I2=imread(‘testpat1.tif’);

subplot(2,2,1);

subimage([0,500],[0,500],RGB1);

subplot(2,2,2);

subimage([0,500],[0,500],RGB2);

subplot(2,2,3);

subimage([0,500],[0,500],I1);

subplot(2,2,4);

subimage([0,500],[0,500],I2);

3．使用不同插值方法放大图像

I=imread(‘rice.tif’);

imshow(I);

I1=imresize(I,1.5,’nearest neighber’);

figure,imshow(I1);

I2=imresize(I,1.5,’bilinear’);

figure,imshow(I2);

I3=imresize(I,1.5,’bicubic’);

figure,imshow(I3);

4．图像的插值旋转

I=imread(‘eight.tif’);

I1=imrotate(I,30,’bilinear’,’crop’);

I2=imrotate(I,30,’nearest’);

imshow(I);

figure,imshow(I1);

figure,imshow(I2);

5．从图像中剪切一块子图

imshow(RGB);

RGB2=imcrop(RGB,[40,50,200,200]);

figure,imshow(RGB2);

（二）设计性实验

读取图像spine.tif，对该图像进行以下处理：

显示原始图像I1；

选用双线性插值方法，将该图像放大2倍，显示放大后的结果I2；

将放大后的图像旋转45度（采用最近邻插值），使旋转后的图像I3与I2大小相同，显示I3；

从图像I3中剪切一块子图，剪切区域为 [60,80,280,320],显示剪切后的图像I4；

I1，I2，I3，I4要求在同一图形窗口显示，并且为每幅图像加上图题。

五、思考题

1. 用montage将两幅图像显示在一齐时，对两幅图像有什么要求？

2． 都有哪些插值方法？各有什么特点？

3． 不同的插值方法对图像放大处理的时间有什么影响？

4． 在选择插值方法时，应考虑哪些因素？

六、实验要求

1、以上实验按顺序进行；

2、设计性实验的实验过程与结果在实验报告中体现出来；

3、实验报告内容还应包括以上4个思考题的回答；

4、实验报告完成。 

实验3 医学图像灰度变换增强

一、实验目的

1、掌握灰度线性变换的原理及实现方法；

2、掌握灰度指数变换的原理及实现方法；

3、掌握灰度对数变换的原理及实现方法；

4、熟悉用MATLAB对医学图像进行灰度变换增强处理。

二、实验原理

当图像成象时曝光不足或过度, 或由于成象设备的非线性和图像记录设备动态范围太窄等因素。都会产生对比度不足的弊病，使图像中的细节分辨不清。这时可将灰度范围线性扩展。

1.线性变换

设f(x,y)灰度范围为[a,b]，g(x,y)灰度范围为[c,d]，那么

（如图1）

2. 分段线性变换

作用：线性扩展感兴趣的灰度范围，抑制不感兴趣的灰度区域。        

设f(x,y)灰度范围为[0,Mf]，g(x,y)灰度范围为[0,Mg],则

（如图2）

图1 灰度线性变换曲线                          图2 灰度分段线性变换曲线

3．指数变换

目的：扩展低灰度区，压缩高灰度区。

变换函数：（a,b,c是参数，可根据需要设置。）

4．对数变换

目的：扩展低灰度区，压缩高灰度区。

变换函数：（a,b,c是参数，可根据需要设置。）

三、实验内容与方法

（一）验证性实验

1、分段线性灰度变换增强

分段线性灰度变换可以突出感兴趣的目标或者灰度区间，相对抑制那些不感兴趣的灰度区域，常用的是三段线性变换。

X1=imread('pout.tif');

figure,imshow(X1);

%定义折点位置

f0=0;g0=0;f1=70;g1=30;f2=180;g2=230;f3=255;g3=255;

%绘制变换曲线

figure,plot([f0,f1,f2,f3],[g0,g1,g2,g3]);

axis tight;xlabel('f');ylabel('g');title('增强变换');

r1=(g1-g0)/(f1-f0);b1=g0-r1*f0;

r2=(g2-g1)/(f2-f1);b2=g1-r2*f1;

r3=(g3-g2)/(f3-f2);b3=g2-r3*f2;

%获取原始图像的长度和宽度

[m,n]=size(X1);X2=double(X1);

%变换矩阵中的每个元素

for i=1:m

    for j=1:n

        %将原始图像矩阵数据赋给f

        f=X2(i,j);

        %初始化变换后的图像矩阵为0

        g(i,j)=0;

        %进行分段灰度变换

if(f>=0)&(f<=f1)

            g(i,j)=r1*f+b1;

elseif(f>=f1)&(f<=f2)

            g(i,j)=r2*f+b2;

elseif(f>=f2)&(f<=f3)

            g(i,j)=r3*f+b3;

        end

    end

end

figure,imshow(mat2gray(g));

2、灰度倒置变换增强

灰度倒置变换就是将原图像的灰度值翻转，黑的变白，白的变黑．

X1=imread('pout.tif');

subplot(1,3,1),imshow(X1);

f1=200;g1=256;

subplot(1,3,2),plot([0,f1],[g1,0]);

axis tight,xlabel('f'),ylabel('g');title('灰度倒置变换');

k=g1/f1;

[m,n]=size(X1);X2=double(X1);

for i=1:m

    for j=1:n

        f=X2(i,j);

        g(i,j)=0;

if(f>=0)&(f<=f1)

            g(i,j)=g1-k*f;

        else

            g(i,j)=0;

        end

    end

end

subplot(1,3,3),imshow(mat2gray(g));

3、灰度对数变换增强

灰度对数变换可以对低灰度区进行较大的扩展而对高灰度区进行压缩，它能使图像灰度的布均匀，与人眼的视觉特性相匹配。

X1=imread('circuit.tif');

subplot(1,3,1),imshow(X1);

c=255/log(256);x=0:1:255;y=c*log(1+x);

subplot(1,3,2),plot(x,y);

axis tight,xlabel('f'),ylabel('g');title('对数变换');

[m,n]=size(X1);

X2=double(X1);

for i=1:m

    for j=1:n

        g(i,j)=c*log(X2(i,j)+1);

    end

end

subplot(1,3,3),imshow(mat2gray(g));

（二）设计性实验

1．读取一副医学图像，用Matlab编程实现其4段灰度变换增强，合理设置变换曲线，使其视觉效果达到最优化。要求在同一图形窗口中同时显示：原始图像、原始图像直方图、变换增强后的图像、变换增强后图像的直方图；在另一图形窗口显示变换曲线。

2．用MATLAB编写程序，实现一幅医学图像的灰度指数变换增强。

四、思考题

1、在灰度分段线性变换中，如何选择各段的折点坐标？

2、灰度倒置是否属于线性变换？

3、指数变换和对数变换有什么异同？

4、指数变换和对数变换分别适用于什么样的图像？

五、实验要求

1、以上实验按顺序进行；

2、实验内容（二）设计性实验的程序与结果须在实验报告中体现出来；

3、实验报告内容还应包括以上4个思考题的回答；

4、实验报告完成，如有30％以上雷同，视为抄袭。 

实验4 医学图像频域滤波增强

一、实验目的

1、掌握频域图像增强的概念、原理与实现过程；

2、掌握频域图像平滑的方法；

3、掌握频域图像锐化的方法。

二、实验原理

频域增强方法将图像从空间域变换到频率域，在图像的频域空间对图像进行滤波处理。傅里叶变换和卷积定理是频域滤波的基础。

在实际应用中，首先需要确定传递函数H(u,v)。如果通过H(u,v)增强F(u,v)的高频信息，如增强图像的边缘信息等，则为高通率波；如果增强F(u,v)的低频信息，如对图像进行平滑操作等，则为低通滤波。频域滤波增强步骤如下：

1）对原始图像f(i,j)进行傅里叶变换得到F(u,v)。

2）将F(u,v)与滤波器函数H(u,v)进行卷及运算得到G(u,v)。

3）对G(u,v)进行傅里叶逆变换，即可求出增强图像g(i,j)。

三、实验内容与方法

1.理想低通滤波器

subplot(1,3,1); imshow(I); J=double(I); f=fft2(J); g=fftshift(f); subplot(1,3,2); imshow(log(abs(g)),[]); [M,N]=size(f); n1=floor(M/2); n2=floor(N/2); d0=5; for i=1:M     for j=1:N         d=sqrt((i-n1)^2+(j-n2)^2); if d<=d0             h=1;         else             h=0;         end         g(i,j)=h*g(i,j);     end end g=ifftshift(g); g=uint8(real(ifft2(g))); subplot(1,3,3); imshow(g);

2.巴特沃斯低通滤波器

subplot(1,3,1); imshow(I); xlabel('a)原始图像'); I2=imnoise(I,'salt & pepper',0.05); subplot(1,3,2); imshow(I2); xlabel('b)含椒盐噪声图像'); f=double(I2); k=fft2(f); g=fftshift(k); [M,N]=size(g); n=2; d0=25; u0=round(M/2); v0=round(N/2); for i=1:M     for j=1:N         d=sqrt((i-u0)^2+(j-v0)^2);         h=1/(1+0.414*d/d0)^(2*n);         y(i,j)=h*g(i,j);        end end y=ifftshift(y); e1=ifft2(y); e2=uint8(real(e1)); subplot(1,3,3);imshow(e2); xlabel('c)巴特沃斯低通滤波器处理结果 ');

3.高斯低通滤波器

subplot(1,3,1); imshow(I); xlabel('a)原始图像'); I2=imnoise(I,'salt & pepper',0.05); subplot(1,3,2); imshow(I2); xlabel('b)含椒盐噪声图像'); f=double(I2); k=fft2(f); g=fftshift(k); [M,N]=size(g); n=2; d0=25; u0=round(M/2); v0=round(N/2); for i=1:M     for j=1:N         d=sqrt((i-u0)^2+(j-v0)^2);         h=exp(-(d/d0)^2);         y(i,j)=h*g(i,j);        end end y=ifftshift(y); e1=ifft2(y); e2=uint8(real(e1)); subplot(1,3,3);imshow(e2); xlabel('c)高通滤波器处理结果');

4.理想高通滤波器

subplot(1,3,1);imshow(I); xlabel('a)原始图像'); J=double(I); f=fft2(J); g=fftshift(f); [M,N]=size(f); n1=floor(M/2); n2=floor(N/2); d0=20; for i=1:M     for j=1:N         d=sqrt((i-n1)^2+(j-n2)^2); if d>=d0             h1=1;             h2=1+0.5;         else             h1=0;             h2=0.5;         end         g1(i,j)=h1*g(i,j);         g2(i,j)=h2*g(i,j);     end end g1=ifftshift(g1); g1=uint8(real(ifft2(g1))); subplot(1,3,2);imshow(g1); xlabel('b)理想高通滤波处理结果'); g2=ifftshift(g2); g2=uint8(real(ifft2(g2))); subplot(1,3,3);imshow(g2); xlabel('b)理想高通加强滤波处理结果');

5.巴特沃斯高通滤波器

subplot(1,3,1);imshow(I); xlabel('a)原始图像'); J=double(I); f=fft2(J); g=fftshift(f); [M,N]=size(f); nn=2; d0=5; m=fix(M/2); n=fix(N/2); for i=1:M     for j=1:N         d=sqrt((i-m)^2+(j-n)^2);         if (d==0)             h=0;         else             h=1/(1+0.414*(d0/d)^(2*nn));         end         s(i,j)=h*g(i,j);     end end s=ifftshift(s); J2=ifft2(s); J3=uint8(real(J2)); subplot(1,2,2);imshow(J3); xlabel('b)巴特沃斯高通滤波处理结果');

6.高斯高通滤波器

subplot(1,2,1);imshow(I); xlabel('a)原始图像'); J=double(I); f=fft2(J); g=fftshift(f); [M,N]=size(f); n=1; d0=10; u0=round(M/2); v0=round(N/2); for i=1:M     for j=1:N         d=sqrt((i-n1)^2+(j-n2)^2);         h=exp(-(d0/d)^n);         y(i,j)=h*g(i,j);             end end y=ifftshift(y); e1=ifft2(y); e2=uint8(real(e1)); subplot(1,2,2);imshow(e2); xlabel('b)高斯高通滤波处理结果');

四、思考题

1、图像增强都有哪些主要方法？哪些属于空域增强法？哪些属于频域增强法？

2、在使用高斯低通滤波器对图像进行平滑处理时，改变滤波器的参数对平滑结果有什么影响？

3、在使用高斯高通滤波器对图像进行锐化处理时，改变滤波器的参数对锐化结果有什么影响？

五、实验与实验报告要求

1、以上实验按顺序进行；

2、在每一实验中，改变滤波器的参数，观察处理结果的变化；

3、实验报告内容还应包括以上3个思考题的回答；

4、实验报告完成，如有30％以上雷同，视为抄袭。 

实验5 医学图像压缩

一、实验目的

1、掌握医学图像压缩编码的原理和基本方法；

2、掌握算术编码的原理及实现方法；

3、掌握行程编码的原理及实现方法。

二、实验原理

算术编码将不同的序列映射到0~1的区域内，该区域表示成可变精度（位数）的二进制小数，越不常见的数据需要的精度越高（更多的位数）。算术编码适合于由相同的重复序列组成的文件，更接近压缩的理论极限，并且不需要码表，只需编码一次就可以得到编码结果。

行程编码将一个相同值的连续串用一个代表值和串长来代替。例如，有一个字符串“aaabccddddd”，经过行程编码后可以“3a1b2c4d”来代替。对图像编码来说，可以定义沿特定方向上具有相同灰度值的相邻像素为一轮，其延续长度称为延续的行程，简称为行程或游程。例如，如果沿水平方向有一串M个像素具有相同的灰度N，则行程编码后，只传递两个值（N,M）就可以代替M个像素的M个灰度值N。

三、实验内容

1、算术编码

新建一个名称为ssbm.m的Matlab脚本文件。输入以下代码：

function ssbm(str) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %算术编码过程 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% l=0;r=1;d=1; %初始间隔 about={ ...     '本实例说明：============'     '字符串不能太长，程序不加判断，请注意溢出；'     '本示例只限定少数字符串；'     '本示例只说明一下算术编码的过程。'} disp(about); %程序限定字符为：a,b,c,d,e p=[0.2 0.3 0.1 0.15 0.25]; %字符的概率分布，sum(p)=1 n=length(str); disp('a b c d e'); disp(num2str(p)); for i=1:n     switch str(i)         case 'a'             m=1;         case 'b'             m=2;         case 'c'             m=3;         case 'd'             m=4;         case 'e'             m=5;         otherwise             error('请不要输入其他字符!');     end     %判断字符     pl=0;  pr=0;     for j=1:m-1         pl=pl+p(j);     end     for j=1:m         pr=pr+p(j);     end     %概率统计     l=l+d*pl;     r=l+d*(pr-pl);     strl=strcat('输入第',int2str(i),'符号的间隔左边界：');     disp(strl);     format long     disp(l);     strr=strcat('输入第',int2str(i),'符号的间隔右边界：');     disp(strr);     disp(r);     d=r-l; end

在Matlab命令行窗口输入：

>>ssbm(‘aedbeecead’)

观察运行结果。改变输入参数，观察并对比运行结果的变化。

2、行程编码

新建一个名称为runlength.m的Matlab脚本文件。输入以下代码并保存：

[zip,info]=xcbm(I);%编码 unzip=dxcbm(zip,info);%解码 %显示原始图像和经编码解码后的图像，显示压缩比 subplot(1,2,1);imshow(I); xlabel('原始图像'); subplot(1,2,2);imshow(unzip); xlabel('行程编码'); cr=info.ratio whos I unzip zip

新建一个名称为xcbm.m的Matlab脚本文件。输入以下代码并保存：

function [zip,info]=xcbm(vector) [m,n]=size(vector); vector=vector(:)'; vector=uint8(vector(:)); L=length(vector); c=vector(1); e(1,1)=c; e(1,2)=0; t1=1; for j=1:L     if((vector(j)==c))         e(t1,2)=double(e(t1,2))+1;     else         c=vector(j);         t1=t1+1;         e(t1,1)=c;         e(t1,2)=1;     end end zip=e; info.rows=m; info.cols=n; [m,n]=size(e); info.ratio=m*n/(info.rows*info.cols);

新建一个名称为dxcbm.m的Matlab脚本文件。输入以下代码并保存：

function unzip=dxcbm(zip,info) zip=uint8(zip); [m,n]=size(zip); unzip=[]; for i=1:m     section=repmat(zip(i,1),1,double(zip(i,2)));     unzip=[unzip section]; end unzip=reshape(unzip,info.rows,info.cols);

运行runlength.m，观察压缩前后图像的变化。

四、思考题

1、算术编码和行程编码各有什么特点？

2、进行实验内容2时，选择5~10幅不同的图像进行实验，观察压缩比的变换情况。总结行程编码的适用对象。

3、算术编码是有损压缩还是无损压缩？行程编码呢？

五、实验与实验报告要求

1、以上实验按顺序进行；

2、在每一实验中，试着改变程序中的参数和输入图像，观察处理结果的变化；

3、实验报告内容还应包括以上3个思考题的回答；

4、实验报告完成，如有30％以上雷同，视为抄袭。 

实验6 医学图像分割

一、实验目的

1、掌握阈值分割的原理与实现方法；

2、掌握微分算子边缘检测的原理及实现方法；

3、掌握各种边缘检测算子的特点及适用范围。

二、实验原理

阈值分割是一种简单有效的图像分割方法。它对物体与背景有较强对比的图像分割特别有效，所有灰度值大于或等于阈值的像素被判断成属于物体，灰度值用255表示前景；否则，这些像素点被排除在物体区域之外，灰度值用0表示背景。

边缘检测是图像分割、目标识别、区域形状提取等图像分析领域十分重要的基础，也是图像识别中提取图像特征的一个重要属性。利用边缘检测来分割图像，其基本思想就是先检测出图像中的边缘点，再按照某种策略将边缘点连接成轮廓，从而构成分割区域。

三、实验内容

1、阈值分割

新建一个名称为threshold.m的Matlab脚本文件。输入以下代码并保存：

imhist(I); %显示原始图像的直方图 %根据上面的直方图选择阈值150，划分图像的前景和背景 newI=im2bw(I,150/255); figure; subplot(1,2,1);imshow(I); xlabel('原始图像'); subplot(1,2,2);imshow(newI); xlabel('阈值法分割图像');

运行threshold.m。观察图像的直方图，判断阈值的选择是否合理。选择不同的图像重复该过程，总结阈值分割法的适用条件。

2、Roberts算子边缘检测

新建一个名称为roberts.m的Matlab脚本文件。输入以下代码并保存：

subplot(1,3,1);imshow(I); xlabel('原始图像'); %以自动阈值选择法对图像进行Roberts算子检测 [BW1,thresh1]=edge(I,'roberts'); %返回当前Roberts算子边缘检测的阈值 disp('Roberts算子自动选择的阈值为') disp(thresh1) subplot(1,3,2);imshow(BW1); xlabel('自动阈值的Roberts算子边缘检测'); %以阈值0.06对图像进行Roberts算子检测 BW2=edge(I,'roberts',0.06); subplot(1,3,3);imshow(BW2); xlabel('阈值为0.06的Roberts算子边缘检测');

运行roberts.m。

3、Sobel算子边缘检测

新建一个名称为sobel.m的Matlab脚本文件。输入以下代码并保存：

subplot(2,2,1);imshow(I); xlabel('原始图像'); %以自动阈值选择法对图像进行Sobel算子检测 [BW1,thresh]=edge(I,'sobel'); %返回当前Sobel算子边缘检测的阈值 disp('Sobel算子自动选择的阈值为:') disp(thresh) subplot(2,2,2);imshow(BW1); xlabel('自动阈值的Sobel算子边缘检测'); %以阈值0.06水平方向对图像进行Sobel算子检测 BW1=edge(I,'sobel',0.06,'horizontal'); subplot(2,2,3);imshow(BW1); xlabel('阈值为0.06水平方向的Sobel算子'); %以阈值0.06垂直方向对图像进行Sobel算子检测 BW2=edge(I,'sobel',0.06,'vertical'); subplot(2,2,4);imshow(BW2); xlabel('阈值为0.06垂直方向的Sobel算子');

运行sobel.m。

4、Canny算子边缘检测

新建一个名称为canny.m的Matlab脚本文件。输入以下代码并保存：

subplot(1,3,1);imshow(I); xlabel('原始图像'); %以自动阈值选择法对图像进行Canny算子检测 [BW,thresh]=edge(I,'canny'); %返回当前Canny算子边缘检测的阈值 disp('Canny算子自动选择的阈值为') disp(thresh) subplot(1,3,2);imshow(BW); xlabel('自动阈值的Canny算子边缘检测'); %以阈值[0.1 0.5]对图像进行Canny算子检测 BW1=edge(I,'canny',[0.1 0.5]); subplot(1,3,3);imshow(BW1); xlabel('阈值为[0.1 0.5]的Canny算子边缘检测');

运行canny.m。

四、思考题

1、阈值分割适用于什么样的图像？

2、什么是边缘？图像边缘有哪些类型？

3、常用的边缘检测算子有哪些？各有什么优缺点？

五、实验与实验报告要求

1、以上实验按顺序进行；

2、在每一实验中，试着改变程序中的参数和输入图像，观察处理结果的变化；

3、实验报告内容还应包括以上3个思考题的回答；

4、实验报告完成，如有30％以上雷同，视为抄袭。