湖南大学医学影像学实验之matlab图像处理实验报告

实验一、空域图像处理

1、灰度线性变换：

I=imread('trees.tif') ;

figure(1) ； imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

J=3*J+74 ; %对J进行线性变换

J=uint8(J) ；

figure(2) ； imshow(J) ; title('线性变换') ;

图像：

思考题：设定不同的斜率值和截距，显示效果会怎样？

答：斜率增加，像素点的灰度值会根据本身的灰度值按比例增加，所以深色部分会变少，浅色部分会变多，整张图片上白色区域会变多；截距增加，所有像素点的灰度值都会增加，所以整张图片会变淡。

添加噪声：

I=imread('trees.tif') ;

M=imnoise(I,'salt & pepper',0.02) ; %添加椒盐噪声

%因为每次添加椒盐噪声都是这两句代码，所以就不重复打出来了，直接看效果图下同～(￣▽￣～)~~~

图像：

椒盐噪声+线性变换：

椒盐噪声是随机产生的噪声，包括高灰度和低灰度的噪声。

线性变换基本上没有多少去噪的功能，更多的应该是用于改变对比度。

(中值)直方图均衡化：

I=imread('rice.png') ;

subplot(2,2,1) ; imshow(I) ; title('原图') ;

subplot(2,2,3) ; imhist(I) ;

J=double(I) ; %把I变成双精度并赋值给J

max=J(1,1) ; min=J(1,1) %把J(1,1)赋值给max和min

[N1,N2]=size(I) ; %得到矩阵I的行和列

for i=1:N1

for j=1:N2

if J(i,j)>max

max=J(i,j) ; %如果元素值大于max，则把元素值赋给max

end

if J(i,j)

min=J(i,j) ; %如果元素值小于min，则把元素值赋给min

end

end

end

n=round((max+min))/2 ; %取元素最大值和最小值的中间值，即中间灰度值

a=(255-max)/(max-n) ；

b=min/(n-min) ；

%根据原图中max和min，求得运算倍率a和b，以确保在接下来的运算中，原图像的max和min可以准确被定为到255和0

for i=1:N1

for j=1:N2

if J(i,j)>=n

J(i,j)=J(i,j)+a*(J(i,j)-n) ;

end

%当元素灰度值大于或等于中间灰度值时，将该元素的灰度值变大

if J(i,j)

J(i,j)=J(i,j)-b*(n-J(i,j)) ;

end

%当元素灰度值小于中间灰度值时，将该元素的灰度值变小

end

end

K=uint8(real(J)) ;

subplot(2,2,2) ; imshow(K) ;title('直方图均衡化') ;

subplot(2,2,4) ; imhist(K) ;

思考题：直方图均衡化是什么意思？它的主要用途是什么？

答：直方图均衡化是将一副像素灰度级范围较窄的图像的像素灰度级的范围扩大并分布均匀。即把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。这样就扩展像原取值的动态范围，提高了对比度和灰度色调的变化，使图像更加清晰。

图像：

1、（中值）直方图均衡化：

2、椒盐噪声+（中值）直方图均衡化：

在添加椒盐噪声后再进行直方图均衡化的过程中，要将添加噪声的过程放在取得运算倍率a和b之后。因为椒盐噪声本身是随机灰度值，如果在计算出之前就添加噪声的话，噪声的灰度值也会被计算在内，会影响原图的运算倍率a和b，这样就会影响直方图均衡化的效果；如果存在灰度值为255和0的噪声，那么直方图均衡化的效果就完全被消除了。所以一定要在之后！！！

均值滤波：

I=imread('rice.png') ;

subplot(1,2,1) ; imshow(I) ; title('原图') ;

[N1, N2]=size(I) ; %得到矩阵I的行和列

I1=double(I) ; %把I变成双精度并赋值给I1

I2=I1 ; %把I1的值都赋值到I2中，即I2是I1的复制

for i=2:N1-1

for j=2:N2-1

A=I1(i-1:i+1,j-1:j+1) ;

%对待处理的像素确定一个模板，该模板包括了其周围的临近像素

s=(A(1,1)+A(1,2)+A(1,3)+A(2,1)+A(2,2)+A(2,3)+A(3,1)+A(3,2)+A(3,3))/9 ;

%得到模板中的全体像素的均值

I2(i,j)=s ; %将所得的均值去替代I2中的同样的位置

end

end

d=uint8(I2) ;

subplot(1,2,2) ; imshow(d) ; title('均值滤波') ;

思考题：

均值滤波的模板大小对处理效果有什么影响？

答：模板越大，可以被参考的数据就越多，所得到的平均值就越真实，去噪效果更好。

图像：

1、均值滤波：

2、添加噪声+均值滤波:

由图可以看出，在没有噪声存在时，均值滤波的效果也能很容易被我们所观察到，不过缺点是，图像本身的细节被模糊了。然后，均值滤波对椒盐噪声的确有去噪的效果，但效果不是和好，很多地方还残存着很多较弱的噪声点。

中值滤波：

I=imread('trees.tif') ;

subplot(1,2,1) ; imshow(I) ; title('原图') ;

[N1, N2]=size(I) ; %得到矩阵I的行和列

I1=double(I) ; %把I变成双精度并赋值给I1

I2=I1 ; %把I1的值都赋值到I2中，即I2是I1的复制矩阵

for i=2:N1-1

for j=2:N2-1

A1=I1(i-1,j-1:j+1) ;

A2=I1(i,j-1:j+1) ;

A3=I1(i+1,j-1:j+1) ;

A=[A1 A2 A3] ;

%对待处理的像素确定一个模板，该模板包括了其周围的临近像素

for p=1:8

for q=1:8

if A(q)

a=A(q) ;

A(q)=A(q+1) ;

A(q+1)=a ;

end

%排序，最终得到从大到小的排序

end

end

I2(i,j)=A(5) ; %用中间值去替代I2中的同样的位置

end

end

d=uint8(I2) ;

subplot(1,2,2) ; imshow(d) ; title('中值滤波') ;

均值滤波，中值滤波的模板大小对处理效果有什么影响？

答：模板越大，可以被参考的数据越多，可以选择在更大范围中的中间值，使单独存在的噪声点更容易被周围的中间值给替代，去噪效果更好。

图像：

1、中值滤波:

细节图：

(与左图相比，可见右图像素点较均匀，不粗糙)

2、添加噪声+中值滤波:

由图可见，在没有噪声存在时，中值滤波的效果不明显，我们很难直接用肉眼观察，只能通过放大了的细节图才能看到效果。而在添加椒盐噪声后，我们会发现去噪效果很明显，椒盐噪声基本上被完全去除了。与均值滤波相比较，中值滤波更适合去除椒盐噪声，效果显著。

拉普拉斯锐化：

I=imread('cameraman.tif') ;

subplot(1,2,1) ; imshow(I) ; title('原图') ;

[N1, N2]=size(I) ; %得到矩阵I的行和列

I1=double(I) ; %把I变成双精度并赋值给I1

I2=I1 ; %把I1的值都赋值到I2中，即I2是I1的复制矩阵for i=2:N1-1

for j=2:N2-1

A=I1(i-1:i+1,j-1:j+1) ;

%对待处理的像素确定一个模板，该模板包括了其周围的临近像素

B=[-1 -1 -1;-1 8 -1;-1 -1 -1] ; %拉普拉斯算子

C=A.*B ;

a=abs(C(1,1)+C(1,2)+C(1,3)+C(2,1)+C(2,2)+C(2,3)+C(3,1)+C(3,2)+C(3,3)) ; %拉普拉斯算子首先将自身与周围的8个像素相减，表示自身与周围像素的差异 I2(i,j)=A(5)-a ;

%将差异加上自身作为新像素的灰度，用它去替代I2中的同样的位置

if a<200

I2(i,j)=255 ;

end

%将差异过小的置为白色，方便对比

end

end

J=uint8(I2) ;

subplot(1,2,2) ; imshow(J) ; title('拉普拉斯') ;

图像：

1、拉普拉斯锐化：

2、添加噪声+拉普拉斯锐化：

由图可以看出，拉普拉斯通过与周围像素点灰度值的差异来强调图像的边缘的细节。但是在去噪方面，拉普拉斯锐化应该说是帮了倒忙。因为噪声本身就是随机产生的像素点，与周围像素点无任何连续性的关系，差异一般都是较大的。所以拉普拉斯锐化在突出边缘和细节的过程中，也顺便增强了噪声。

实验二、傅立叶变换

高斯滤波：

I=imread('cameraman.tif') ;

fftI=fft2(I) ; %二维离散傅里叶变换

sfftI=fftshift(fftI) ; %直流分量移到频谱中心

RR=real(sfftI) ; %取傅里叶变换的实部

II=imag(sfftI) ; %取傅里叶变换的虚部

A=sqrt(RR.^2+II.^2) ; %计算频谱幅值

A=(A-min(min(A)))/(max(max(A))-min(min(A)))*255 ; %归一化

subplot(1,2,1) ; imshow(A) ; title('傅里叶变换') ;

J=fspecial('gaussian') ; %调用高斯函数

K=imfilter(A,J) ; %用高斯函数对A进行滤波

subplot(1,2,2) ; imshow(K) ; title('高斯滤波') ;

图像：

细节图：（与左图相比，可见右图中心些许像素点已被滤去）

思考题：

1、傅里叶变换有哪些重要的性质?

答：线性，对偶性，平移性，尺度变换，微分关系，时域、空域卷积定理。2、图像的二维频谱在显示和处理时应注意什么?

答：1、进行傅里叶变换的图像应是灰度图形，rgb彩色图像无法进行变换；

2、使用fftshift函数将频谱的零频分量移至频谱的中心；

3、要进行归一化。

I1=imread('rice.png') ;

subplot(1,2,1) ; imshow(I1) ; title('原图') ;

m=double(I1) ; %把I变成双精度并赋值给m

f=fft2(m) ; %二维离散傅里叶变换

f=fftshift(f) ; %直流分量移到频谱中心

[N1,N2]=size(f) ; %得到矩阵f的行和列

n1=round(N1/2) ; %取行的中点并取整

n2=round(N2/2) ; %取列的中点并取整

d0=50 ; %设定临界值

for i=1:N1

for j=1:N2

d=sqrt((i-n1)^2+(j-n2)^2) ; %计算频谱幅值

if d<=d0

h=0.5 ; %低于或等于临界值，本应滤去，即设定倍率h=0

但为了提高对比度，设定h=0.5

else h=1 ; %高于临界值，保留，即设定倍率h=1

end

y(i,j)=h*f(i,j); %将原灰度值乘以倍率得到新的灰度值

end

end

y=ifftshift(y) ; %反中心平移

A=ifft2(y) ; %二维傅里叶反变换

B=uint8(real(A)) ;

subplot(1,2,2) ; imshow(B) ; title('高通滤波') ;

图像：

由图可见，高通滤波滤去了低频信号，使图像变暗了很多，明显的突出了图像的边缘和细节，这是图像锐化的一种方式。

I1=imread('rice.png') ;

subplot(1,2,1) ; imshow(I1) ; title('原图') ;

m=double(I1) ; %把I变成双精度并赋值给m

f=fft2(m) ; %二维离散傅里叶变换

f=fftshift(f) ; %直流分量移到频谱中心

[N1,N2]=size(f) ; %得到矩阵f的行和列

n1=round(N1/2) ; %取行的中点并取整

n2=round(N2/2) ; %取列的中点并取整

d0=50 ; %设定临界值

for i=1:N1

for j=1:N2

d=sqrt((i-n1)^2+(j-n2)^2) ; %计算频谱幅值

if d<=d0

h=1 ; %高于临界值，保留，即设定倍率h=1

else h=0 ; %低于或等于临界值，滤去，即设定倍率h=0

end

y(i,j)=h*f(i,j); %将原灰度值乘以倍率得到新的灰度值

end

end

y=ifftshift(y) ; %反中心平移

A=ifft2(y) ; %二维傅里叶反变换

B=uint8(real(A)) ;

subplot(1,2,2) ; imshow(B) ; title('高通滤波') ;

图像：

有图可见，低通滤波滤去了高频信号，使得图像的边缘都被弱化了，图像整体变模糊了。这是图像平滑的一种方式。

实验三、图像的几何变换

平移：

I=imread('trees.tif') ;

figure(1); imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

[N1,N2]=size(J) ; %得到矩阵J的行和列

K=zeros(N1,N2) ; %设定一个N1XN2的零矩阵K

A=[1 0 60;0 1 30;0 0 1] ; %设定一个平移矩阵A

for i=1:N1

for j=1:N2 %预计平移后超出原矩阵范围的数据将不予传输

B=[i;j;1] ; %将初始坐标设置为矩阵B

C=A*B ; %经过矩阵运算得到矩阵C

x=C(1) ; y=C(2) ; %把C中的元素值、即平移后的坐标赋给K

K(x,y)=J(i,j) ; %将J矩阵中该点的值赋给K矩阵中平移后的对应点

end

end

figure(2) ; imshow(K) ; title('平移') ;

图像：

由图可见，平移后的图像，超出原矩阵范围的数据会丢失，所以在编程过程中，超过范围的数据就不予传输。如果没有丢弃那些数据，就会导致平移后的图像会比原图像大。

思考题：改变水平和垂直的偏移量，观察显示？

答：水平方向和垂直方向的平移距离会改变，图像中黑色区域会改变，图像数据丢失的量也会改变。

I=imread('trees.tif') ;

subplot(1,2,1) ; imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

[N1,N2]=size(J) ; %得到矩阵J的行和列

K=zeros(N1,N2) ; %设定一个N1XN2的零矩阵K

A=[1 0 0;0 -1 0;0 0 1] ; %设定一个水平镜像矩阵A

for i=1:N1

for j=1:N2

B=[i;j;1] ; %将初始坐标设置为矩阵B

C=A*B ; %经过矩阵运算得到矩阵C

x=C(1) ;

y=C(2)+N2+1 ; %把C中的元素值、即水平镜像后的坐标赋给K

K(x,y)=J(i,j); %将J矩阵中该点的值赋给K矩阵中水平镜像后的对应点 end

end

K=uint8(K) ;

subplot(1,2,2) ; imshow(K) ; title('水平镜像') ;

图像：

在水平镜像的过程中，起始y坐标经过水平矩阵的运算后得到的最终y坐标会变成负值。而坐标值在检索过程中只能为正整数，所以需要对y坐标进行单独的平移运算。

I=imread('trees.tif') ;

subplot(1,2,1) ; imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

[N1,N2]=size(J) ; %得到矩阵J的行和列

K=zeros(N1,N2) ; %设定一个N1XN2的零矩阵K

A=[-1 0 0;0 1 0;0 0 1] ; %设定一个垂直镜像矩阵A

for i=1:N1

for j=1:N2

B=[i;j;1] ; %将初始坐标设置为矩阵B

C=A*B ; %经过矩阵运算得到矩阵C

x=C(1)+N1+1 ;

y=C(2) ; %把C中的元素值、即垂直镜像后的坐标赋给K

K(x,y)=J(i,j) ; %将J矩阵中该点的值赋给K矩阵中垂直镜像后的对应点 end

end

K=uint8(K) ;

subplot(1,2,2) ; imshow(K) ; title('垂直镜像') ;

图像：

垂直镜像和水平镜像的过程很类似，只不过是起始x坐标经过水平矩阵的运算后得到的最终x坐标会变成负值。所以我们也需要对x坐标进行单独的平移运算。

I=imread('trees.tif') ;

figure(1) ; imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

[N1,N2]=size(J) ; %得到矩阵J的行和列

K=zeros(N1/2,N2/2) ; %设定一个(N1/2)X(N2/2)的零矩阵K

for i=2:2:N1

for j=2:2:N2

x=i/2 ; y=j/2 ; %将J矩阵的坐标值减半并赋给K的坐标值

K(x,y)=J(i,j) ; %将J矩阵中该点的值赋给K矩阵中缩小后的对应点

end

end

K=uint8(K) ;

figure(2) ; imshow(K) ; title('缩小') ;

图像：

思考题：改变缩小比例，看看效果如何？

答：缩小，本质是根据一定的比例去丢失一部分数据，剩下一部分数据。剩下的数据维持缩小之前的大概图像；如果缩小的越小，那么需要丢失的数据也越多，剩下的数据也就会越少，这样一来图像也会越模糊。

而且因为坐标的原因、我编写的程序让图像缩小和放大的倍数只能是整倍数。

I=imread('trees.tif') ;

figure(1) ; imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

[N1,N2]=size(J) ; %得到矩阵J的行和列

K=zeros(2*N1,2*N2) ; %设定一个(2*N1)X(2*N2)的零矩阵K

A=[2 0 0;0 2 0;0 0 1] ; %设定一个放大矩阵A

for i=1:N1

for j=1:N2

B=[i;j;1] ; %将初始坐标设置为矩阵B

C=A*B ; %经过矩阵运算得到矩阵C

x=C(1) ; y=C(2) ; %把C中的元素值、即放大后的坐标赋给K

K(x,y)=J(i,j) ; %将J矩阵中该点的值赋给K矩阵中放大后的对应点

end

end

for i=1:2*N1

for j=1:2*N2

if K(i,j)==0

K(i,j)=30*rand(1) ;

end

%对于未填充的数值都是0，为了填充数值，用随机数填充

end

end

K=uint8(K) ;

figure(2) ; imshow(K) ; title('放大') ;

图像：

思考题：改变缩放比例，看看效果如何？

答：放大，本质是除了本身的数据外，还根据一定的规律去填充数据；放大的越大，填充值所占的比例会越大，本身的数据所占比例就越小，图像也就会越模糊。

I=imread('trees.tif') ;

figure(1); imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

[N1,N2]=size(J) ; %得到矩阵J的行和列

K=zeros(N1,N2) ; %设定一个(N1/2)X(N2/2)的零矩阵K

b=pi/4 ; %设定一个旋转角度

A=[cos(b) sin(b) 0;-sin(b) cos(b) 0;0 0 1] ; %设定一个旋转矩阵A

for i=1:N1

for j=1:N2

B=[i;j;1] ; %将初始坐标设置为矩阵B

C=A*B ; %经过矩阵运算得到矩阵C

x=round(C(1)) ;

y=round(C(2))+200 ; %把C中的元素值、即旋转后的坐标赋给K

K(x,y)=J(i,j) ; %将J矩阵中该点的值赋给K矩阵中旋转后的对应点

end

end

K=uint8(K) ;

figure(2) ; imshow(K) ; title('旋转') ;

图像：

坐标经过旋转矩阵运算后得到的值包含非整数，但是由于是坐标的原因，可以被检索的只能是整倍数。所以我的处理方式是对坐标进行平方运算，这样就可以使坐标被检索了；但是，缺点是原图会被放大，所以还不是完美的旋转。

思考题：改变旋转角度，显示效果会怎么样？

答：旋转，是图像以原图像中心点为原点旋转。改变旋转角度后，图像的旋转角度会改变，旋转后的整张图片的大小也会有变化，但是图片中的原图大小不变。

实验四、图像分割

二值化：

I=imread('trees.tif') ;

subplot(1,2,1); imshow(I) ; title('原图') ;

J=double(I) ; %把I变成双精度并赋值给J

max=J(1,1) ; min=J(1,1) ; %把J(1,1)赋值给max和min

[N1,N2]=size(I) ; %得到矩阵I的行和列

for i=1:N1

for j=1:N2

if J(i,j)>max

max=J(i,j) ; %如果元素值大于max，则把元素值赋给max

end

if J(i,j)

min=J(i,j) ; %如果元素值小于min，则把元素值赋给min

end

end

end

n=round((max+min))/2 ; %取元素最大值和最小值的中间值，即中间灰度值为阈值

K=zeros(N1,N2) ; %设定一个N1XN2的零矩阵K

for i=1:N1

for j=1:N2

if J(i,j)<=n

K(i,j)=0 ;

else

K(i,j)=1 ;

end

%利用n这个阈值，将灰度值小于或的等于它的都设定为黑，大于它的设定为白 end

end

subplot(1,2,2) ; imshow(K) ; title('二值化') ;

图像：