医学成像与图像处理实验报告.docx
《医学成像与图像处理实验报告.docx》由会员分享，可在线阅读，更多相关《医学成像与图像处理实验报告.docx（17页珍藏版）》请在冰豆网上搜索。
医学成像与图像处理实验报告
医学成像与图像处理实验报告
实验一空间域图像增强
图像：
Eimage-007.img
1.1平滑处理：
分别用5x5,7x7的平滑模板作平滑处理。
clear;
fid=fopen（'Eimage-007.img','rb'）;
C=fread（fid,65536,'float64'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
%A=A';
h1=fspecial（'average',5）;%5x5平滑模板
h2=fspecial（'average',7）;%7x7平滑模板
A1=imfilter（A,h1）;
A2=imfilter（A,h2）;
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
subplot（221）;imshow（A,[minmin,maxmax]）;title（'原始图像'）;
maxmax=max（max（A1））;
minmin=min（min（A1））;
subplot（223）;imshow（A1,[minmin,maxmax]）;title（'5x5平滑处理图像'）;
maxmax=max（max（A2））;
minmin=min（min（A2））;
subplot（224）;imshow（A2,[minmin,maxmax]）;title（'7x7平滑处理图像'）;
实验结果：
原始图像经过平滑处理后变得模糊，并且7x7平滑处理的图像比5x5平滑处理的图像更模糊。
这说明领域半径越大，模糊程度就越大。
1.2用罗伯特Robert梯度法提取图像边缘。
clear;
fid=fopen（'Eimage-007.img','rb'）;
C=fread（fid,65536,'float64'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
%A=A'
grayPic=mat2gray（A）;%图像矩阵的归一化操作
[m,n]=size（grayPic）;
newGrayPic=grayPic;%为保留图像边缘的一个像素
robertsNum=0;%经roberts算子计算得到的每个像素的值
robertThreshold=0.1;%设定阈值
forj=1:
m-1%进行边界提取
fork=1:
n-1
robertsNum=abs（grayPic（j,k）-grayPic（j+1,k+1））+...
abs（grayPic（j+1,k）-grayPic（j,k+1））;
if（robertsNum>robertThreshold）
newGrayPic（j,k）=255;
else
newGrayPic（j,k）=0;
end
end
end
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
subplot（121）;imshow（A,[minmin,maxmax]）;title（'原始图像'）;
maxmax=max（max（newGrayPic））;
minmin=min（min（newGrayPic））;
subplot（122）;imshow（newGrayPic,[minmin,maxmax]）;title（'Robert算子处理后图像'）;
实验结果：
罗伯特Robert梯度法边缘检测，Roberts算子采用对角线方向相邻两像素之差近似的梯度幅值来检测边缘。
该算子定位较准确，但对噪声比较敏感，检测水平和竖直边缘效果好于斜向边缘。
1.3用拉普拉斯算子做锐化处理
clear;
fid=fopen（'Eimage-007.img','rb'）;
C=fread（fid,65536,'float64'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
subplot（221）;imshow（A,[minmin,maxmax]）;title（'原始图像'）;
fGray=A;
[r,c]=size（fGray）;
fLapEdge=double（fGray）;
M=fLapEdge;
fori=2:
r-1
forj=2:
c-1
M（i,j）=4*fGray（i,j）-fGray（i-1,j）-fGray（i,j+1）-...
fGray（i+1,j）-fGray（i,j-1）;%四领域模板
end
end
s=M+fGray;
S=uint8（s）;
maxmax=max（max（M））;
minmin=min（min（M））;
subplot（222）;imshow（M,[minmin,maxmax]）;title（'拉普拉斯算子处理后的图像'）;
maxmax=max（max（S））;
minmin=min（min（S））;
subplot（223）;imshow（S,[minmin,maxmax]）;title（'锐化后的图像'）;
实验结果：
用拉普拉斯算子四领域模板对图像进行了锐化处理，增大领域间像素的差值，加强了图像的边缘和轮廓。
实验二频率域图像增强
图像：
Eimage-007.img
2.1频率域平滑处理：
理想低通和巴特沃斯低通。
clear;
fid=fopen（'Eimage-007.img','rb'）;
C=fread（fid,65536,'float64'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
subplot（221）;imshow（A,[minmin,maxmax]）;title（'原始图像'）;
Fimg=fft2（A）;%傅里叶变换
Fimg=fftshift（Fimg）;%将变换的原点移到频率矩形的中心
[M,N]=size（Fimg）;
n1=floor（M/2）;%取整
n2=floor（N/2）;%取整
%半径为5的理想低通滤波处理
dist1=5;
fori=1:
M
forj=i:
N
if（sqrt（（（i-n1）^2+（j-n2）^2））<=dist1）
h=1;
else
h=0;
end
Fimg（i,j）=h.*Fimg（i,j）;
end
end
g1=ifftshift（Fimg）;%反移中
img1=real（ifft2（g1））;
maxmax=max（max（img1））;
minmin=min（min（img1））;
subplot（223）;imshow（img1,[minmin,maxmax]）;title（'半径5的理想低通滤波'）;
%巴特沃斯低通滤波
Fimg=fft2（A）;%傅里叶变换
Fimg=fftshift（Fimg）;%将变换的原点移到频率矩形的中心
[M,N]=size（Fimg）;
n1=floor（M/2）;%取整
n2=floor（N/2）;%取整
n=1;%对n赋初值
%半径为5的巴特沃斯低通滤波处理
d0=5;
fori=1:
M
forj=i:
N
d=sqrt（（（i-n1）^2+（j-n2）^2））;%到傅里叶变换中心的距离
h=1/（1+（d/d0）^（2*n））;%滤波函数
Fimg（i,j）=h.*Fimg（i,j）;
end
end
g1=ifftshift（Fimg）;%反移中
img1=real（ifft2（g1））;
maxmax=max（max（img1））;
minmin=min（min（img1））;
subplot（224）;imshow（img1,[minmin,maxmax]）;title（'半径5的巴特沃斯低通滤波'）;
实验结果：
分别用理想低通滤波和巴特沃斯低通滤波进行频率域平滑处理，由于理想低通滤波器有陡峭频率的截止特性，会产生明显的振铃现象使图像变得模糊，巴特沃斯的振铃现象不明显。
2.2频率域锐化处理：
理想高通和巴特沃斯高通。
clear;
fid=fopen（'Eimage-007.img','rb'）;
C=fread（fid,65536,'float64'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
subplot（221）;imshow（A,[minmin,maxmax]）;title（'原始图像'）;
Fimg=fft2（A）;%傅里叶变换
Fimg=fftshift（Fimg）;%将变换的原点移到频率矩形的中心
[M,N]=size（Fimg）;
n1=floor（M/2）;%取整
n2=floor（N/2）;%取整
%半径为5的理想高通滤波处理
dist1=5;
fori=1:
M
forj=i:
N
if（sqrt（（（i-n1）^2+（j-n2）^2））>=dist1）
h=1;
else
h=0;
end
Fimg（i,j）=h.*Fimg（i,j）;
end
end
g1=ifftshift（Fimg）;%反移中
img1=real（ifft2（g1））;
maxmax=max（max（img1））;
minmin=min（min（img1））;
subplot（223）;imshow（img1,[minmin,maxmax]）;title（'半径5的理想高通滤波'）;
Fimg=fft2（A）;%傅里叶变换
Fimg=fftshift（Fimg）;%将变换的原点移到频率矩形的中心
[M,N]=size（Fimg）;
n1=floor（M/2）;%取整
n2=floor（N/2）;%取整
n=1;%对n赋初值
%半径为5的巴特沃斯高通滤波处理
d0=5;
fori=1:
M
forj=i:
N
d=sqrt（（（i-n1）^2+（j-n2）^2））;%到傅里叶变换中心的距离
h=1/（1+（d0/d）^（2*n））;%滤波函数
Fimg（i,j）=h.*Fimg（i,j）;
end
end
g1=ifftshift（Fimg）;%反移中
img1=real（ifft2（g1））;
maxmax=max（max（img1））;
minmin=min（min（img1））;
subplot（224）;imshow（img1,[minmin,maxmax]）;title（'半径5的巴特沃斯高通滤波'）;
实验结果：
采用理想高通和巴特沃斯高通滤波器对图像锐化处理，理想高通有明显振铃现象，图像的边缘有抖动现象；巴特沃斯高通效果较好，H（u,v）是渐变的，振铃现象不明显。
实验三图像恢复
图像：
smooth.img
图像受大气湍流的影响变得模糊，K=0.0025。
用逆滤波的方法恢复原图像（用巴特沃斯滤波器）。
用维纳滤波的方法恢复原图像。
clear;
fid=fopen（'smooth.img','rb'）;
C=fread（fid,65536,'float32'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
A=A';
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
imshow（A,[minmin,maxmax]）;title（'模糊图像'）;
[m,n]=size（A）;
k=0.0025;
foru=1:
m
forv=1:
n
H（u,v）=exp（（-k）*（（（u-m/2）^2+（v-n/2）^2）^（5/6）））;
end
end
F0=fftshift（fft2（A））;
F1=F0./H;
Fimg=fftshift（F1）;%将变换的原点移到频率矩形的中心
[M,N]=size（Fimg）;
n1=floor（M/2）;%取整
n2=floor（N/2）;%取整
p=10;%对p赋初值
%半径为5的巴特沃斯高通滤波处理
d0=5;
fori=1:
M
forj=i:
N
d=sqrt（（（i-n1）^2+（j-n2）^2））;%到傅里叶变换中心的距离
h=1/（1+（d0/d）^（2*p））;%滤波函数
Fimg（i,j）=h.*Fimg（i,j）;
end
end
I2=ifftshift（Fimg）;%反移中
I3=uint8（real（ifft2（I2）））;
figure
（2）
imshow（I3）;title（'逆滤波恢复图像'）;
K=0.1;
foru=1:
m
forv=1:
n
H（u,v）=exp（（-k）*（（（u-m/2）^2+（v-n/2）^2）^（5/6）））;
H0（u,v）=（abs（H（u,v）））^2;
H1（u,v）=H0（u,v）/（H（u,v）*（H0（u,v）+K））;
end
end
F2=H1.*F0;
I4=ifft2（fftshift（F2））;
I5=uint8（I4）;
figure（3）
imshow（I5）;title（'维纳滤波恢复图像'）;
实验结果：
图像复原的过程首先利用退化现象的某种先验知识，建立退化现象的数学模型，然后再根据退化模型进行反向的推演运算，以恢复原来的景物图像。
逆滤波复原方法数学表达式简单，物理意义明确。
但是，当存在噪声时，而且H（u,v）很小或者为零的时候，噪声会被放大，会对复原的图像产生巨大的影响，有可能使恢复的图像和原图相差很大。
维纳滤波的思想是在假想图像信号可以看作平稳随机过程的前提下，按照使恢复的图像与原图像的均方差最小原则来恢复图像。
可以看出维纳滤波在图像受噪声影响时效果比逆滤波要好，而且噪声越强优势越明显。
实验四图像分割
图像：
Eimage-007.img
用区域生长的方法分割出大脑区域。
clear;
fid=fopen（'Eimage-007.img','rb'）;
C=fread（fid,65536,'float64'）;
fclose（fid）;
fori=1:
256
forj=1:
256
A（i,j）=C（（i-1）*256+j）;
end
end
figure
（1）
maxmax=max（max（A））;
minmin=min（min（A））;
imshow（A,[minmin,maxmax]）;title（'原始图像'）;
I=double（A）;%转换为灰度值是0-1的双精度
[M,N]=size（I）;%得到原图像的行列数
[y,x]=getpts;%获得区域生长起始点
x1=round（x）;%横坐标取整
y1=round（y）;%纵坐标取整
seed=I（x1,y1）;%将生长起始点灰度值存入seed中
Y=zeros（M,N）;%作一个全零与原图像等大的图像矩阵Y，作为输出图像矩阵
Y（x1,y1）=1;%将Y中与所取点相对应位置的点设置为白点
sum=seed;%储存符合区域生长条件的点的灰度值的总和
suit=1;%储存符合区域生长条件的点的总个数
count=1;%每次判断一点周围八点符合条件的新点的数目
threshold=40;%域值，即某一点与周围八点的绝对差值要小于阈值
whilecount>0%判断是否有新的符合生长条件的点，若没有，则结束
s=0;%判断一点周围八点时，符合条件的新点的灰度值之和
count=0;
fori=1:
M
forj=1:
N
ifY（i,j）==1
if（i-1）>0&（i+1）<（M+1）&（j-1）>0&（j+1）<（N+1）%判断此点是否为图像边界上的点
foru=-1:
1%判断点周围八点是否符合域值条件
forv=-1:
1%u,v为偏移量
ifY（i+u,j+v）==0&abs（I（i+u,j+v）-seed）<=threshold%判断是否未存在于输出矩阵Y，并且为符合域值条件的点
Y（i+u,j+v）=1;%符合以上两条件即将其在Y中与之位置对应的点设置为白点
count=count+1;%新的、符合生长条件的点的总个数
s=s+I（i+u,j+v）;%新的、符合生长条件的点的总灰度数
end
end
end
end
end
end
end
suit=suit+count;%目前区域所有符合生长条件的点的总个数
sum=sum+s;%目前区域所有符合生长条件的点的总灰度值
seed=sum/suit;%计算新的灰度平均值
end
figure
（2）;imshow（Y）;title（'分割出的大脑区域'）;
实验结果：
区域生长是指从某个像素出发，按照一定的准则，逐步加入邻近像素，当满足一定的条件时，区域生长终止。
区域生长是从某个或者某些像素点出发，最后得到整个区域，进而实现目标的提取。
1、选取图像中的一点为种子点（种子点的选取需要具体情况具体分析）。
2、在种子点处进行8邻域或4邻域扩展，判定准则是：
如果考虑的像素与种子像素灰度值相差的绝对值小于某个门限T，则将该像素包括进种子像素所在的区域。
3、当不再有像素满足加入这个区域的准则时，区域生长停止。
实验五图像重建
图像：
atten.rad
用滤波反投影的方法重建图像。
clear;
fid=fopen（'atten.rad','rb'）;
C=fread（fid,16384,'float32'）;
fclose（fid）;
fori=1:
128
forj=1:
128
A（i,j）=C（（i-1）*128+j）;
end
end
%A=A'
H=zeros（128,128）;
forx1=1:
128
fory1=1:
128
fork=1:
128
b=（-2）/（（（pi*1）^2）*（4*（（y1-k）^2）-1））;%S-L滤波函数公式为-2/（pi2d2（4n2-1）），n为整
H（x1,y1）=A（x1,k）*b+H（x1,y1）;
end
end
end
%反投影
pic=zeros（128,128）;
fori=1:
128
forj=1:
128
x2=i-64.5;%图片中心点为64.5
y2=64.5-j;
forcol=1:
128
angle=（col-1）*2*pi/128;
xr=x2*cos（angle）+y2*sin（angle）;%反投影，公式xr,m=xi*cos（fi）+yj*sin（fi）
xr=xr+64.5;
if（1value_xr=H（col,floor（xr））*（1-（xr-floor（xr）））+H（col,ceil（xr））*（（xr-floor（xr）））;
pic（i,j）=value_xr+pic（i,j）;
end
end
end
end
figure;
maxmax=max（max（A））;
minmin=min（min（A））;
subplot（121）;imshow（A,[minmin,maxmax]）;title（'原始图像'）
maxmax=max（max（pic））;
minmin=min（min（pic））;
subplot（122）;imshow（pic,[minmin,maxmax]）;title（'滤波反投影算法重建后图像'）
实验结果：
用了卷积计算，射束计算与内插，反投影重建的方法实现了滤波反投影重建图像。
使用了S-L滤波函数，S-L滤波器对投影中的高频成分具有抑制作用，进而使重建图像的振荡响应减小，对含噪声的投影数据，重建质量用RL较好，但在低频段效果不及R-L滤波器质量好，重建图像还是存在一定误差。
升级会员 