生物医学工程课程设计报告-MRI图像的有约束恢复和几何失真校正.doc

生物医学工程课程设计报告-MRI图像的有约束恢复和几何失真校正

课程设计报告

题 目：MRI图像的有约束恢复和几何失真校正
专 业： 生物医学工程
班 级
学 号：
姓 名：
指导教师：

天津医科大学 生物医学工程学院
2010年 10月 27日
（一） 开题背景
图像复原就是在研究图象退化原因的基础上，以退化图像为依据，根据一定的先验知识设计一种算子，从而估计出理想场景的操作。一般得到一幅数字化图像后都会先使用图像复原技术进项处理，然后再做增强处理。由于不同应用领域的图像有不同的退化原因，所以对同一幅退化图像，不同应用领域采用不同的复原方法。图像复原可以看成是图像复原的逆过程，对退化复原的一般可采用两种方法：一种方法适用于图像缺乏先验知识的情况下，此时可对退化过程（噪声和模糊）建立模型，进行描述。并寻找一种去除和消弱其影响的过程，从而改善图像质量。另一方面若对于原始图像有足够的先验知识，则对原始图像简历一个数学模型，并根据它对退化图像进行拟合，将其转化为一个检验问题。但大多数情况下退化过程是不可知的，由于图像模糊的同时，噪声和干扰也会同时存在，这也为复原过程带来困难和不确定性。本文针对当前主流的图像复原算法进行分析，并归纳和总结，并进行Matlab的仿真实验，为人们的不同应用场合及不同的图像数据条件下选择不同的复原算法提供了一定依据。

（二）课题目的
1. 提高分析问题和解决问题的能力，进一步巩固数字图像处理系统中的基本原理与方法；
2了解图像退化的原因，并建立图像退化模型；
3.学习用MATLAB图像复原（维纳滤波、约束最小平方滤波）和几何失真校正的方法。

（三）课题研究的主要内容
随着计算机技术和影像诊断技术(MRI、CT、DSA等的发展) ，立体放射治疗技术(X-刀、-刀治疗系统)已广泛运用于颅内小病变的治疗，靶点位置精确度成为治疗的成败关键因素之一，MRI以其优良的软组织对比度、多平面直接成像的优点，使其广泛应用于立体放射治疗靶点定位，但由于MRI成像过程中的图象退化以及几何失真（实际组织空间结构与它在图像上表现存在差别）是影响图像清晰度和靶点位置精确度的主要因素。
在学习了医学图像处理、医学成像系统课程的基础上，利用维纳滤波（最小均方误差）、约束最小平方滤波器、几何失真校正实现MRI图像的恢复。
1. 对给定的MRI图通过模拟水平运动模糊建立退化函数，模糊原图像并加入高斯噪声；
2. 应用维纳（最小均方误差LMS）滤波器恢复图像；
3. 应用约束最小平方滤波器恢复图像；
4. 对图像进行仿射变换，使其变为原来大小的1/4，再分别利用最近邻插值、双线型插值和三次内插法进行集合校正，将其校正为原始大小。
（四）原理和方法
图像在形成、记录、处理和传输过程中，由于成像系统、记录设备、传输介质和处理方法的不完善，从而导致图像质量下降，这种现象就成为图像退化。例如，光学系统中的衍射、光电转换器件的非线性、光学系统中的像差、大气湍流的扰动效应、曝光噪声干扰、图像运动造成的模糊性以及几何畸变等等。
（1） 成像系统的像差、畸变、有限带宽等造成的图像失真；
（2） 涉嫌辐射、大气湍流等造成的照片畸变；
（3） 携带遥感仪器的飞机或卫星运动的不稳定，以及地球自传等因素引起的照片几何失真；
（4） 模拟图像在数字化的过程中，由于会损失掉部分细节，因而造成图像质量下降；
（5） 拍摄时，相机与景物之间的相对运动产生的运动模糊；
（6） 镜头聚焦不准产生的散焦模糊；
（7） 底片感光、图像显示时会造成记录显示失真；
（8） 成像系统中存在的噪声干扰。
图像复原就是要尽可能恢复退化图像的本来面目，它是沿图像退化的逆过程进行处理。
典型的图像复原是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用各种逆退化处理方法进行恢复，使图像质量得到改善。图像复原过程如下：
分析退化原因 建立退化模型 反向推演 反向推演

图像复原的关键在于建立图像退化模型。这个退化模型应该能够反映图像退化的原因。通常将退化原因作为线性系统退化的一个因素来对待，从而建立系统退化模型来近似描述图像函数的退化。如下图所示，这是一种简单的通用退化模型，它将图像的退化过程模型化为一个退化系统（或退化算子）H。由图可见一幅纯净的图像f(x,y)由于通过一个系统H以及引 进外来加性噪声n（x，y）而退化为一幅图像g(x,y)。
n(x,y)
f(x,y) ─→ ─→⊕─→g(x,y)

H

输入输出关系式为：g(x,y)=H[f(x,y)]+n(x,y) (1)
在这里，n(x,y)是一种统计性质的信息。在实际应用中，往往假设噪声是白噪声，即它的频谱密度为常熟，并且与图像不相关。
在对退化系统进行了线性系统和空间不变系统的近似之后，连续函数的退化模型在空域中可以写成：
g(x,y)=f(x,y)*h(x,y)+n(x,y) (2)
在频域中可以写成：
G(u,v)=F(u,v)H(u,v)+N(u,v) (3)
其中，G(u,v)、F(u,v)、N(u,v)分别是退化图像g(x,y)、原图像f(x,y)、噪声信号n(x,y)的傅立叶变换；H(u,v)是系统的点冲击响应函数h(x,y)的傅立叶变换，称为系统在频率域上的传递函数。
可见，图像复原实际上就是已知g(x,y)求f(x,y)的问题或已知G(u,v)求F(u,v)的问题，它们的不同之处在于一个是空域，一个是频域。

图像复原的主要目的是在给定退化图像g(x,y)以及退化函数H和噪声的某种了解或假设时，估计出原始图像f(x,y)。现在的问题是退化函数H一般是不知道的。因此必须在进行图像复原前对退化函数进行估计。

1.几种有效的图像复原方法
（1）维纳滤波复原
维纳滤波（N.Wiener最先在1942年提出的方法）是一种最早，也是最为人们熟知的线性图像复原方法。

f 的相关矩阵

n 的相关矩阵

用块循环矩阵表示

A和B中的元素对应Rf和Rn中的相关元素的傅里叶变换，这些相关元素的傅里叶变换称为图像和噪声的功率谱。
令

是退化图像的傅里叶变换
是退化函数

是退化函数的复共轭
是噪声的功率谱
是原始图像的功率谱
维纳滤波器的传递函数

（2）有约束最小平方滤波方法
有约束最小平方滤波就是利用已知的复原图像和噪声的均值与方差， 通过退化函数，经过反傅立叶变换，求出原始图像。

其中，

为的傅里叶变换。的取值控制着对估计图像所加光滑性约束的强度。P（u,v） 为用 Q（高通卷积滤波运算）实现的高通滤波器的传函，决定了不同频率所受光滑性影响的程度。
2.仿射变换
仿射变换，可以用以下函数来描述：
，其中，A是变形矩阵，b是平移矩阵。
尺度变换：
变换矩阵：,S≥0

3.灰度级插值
实现几何运算有两种方法：
其一为前向映射法，即将输入像素的灰度一个个地转移到输出图像中，如果一个输入像素被映射到四个输出像素之间的位置，则其灰度值就按插值法在四个输出像素之间进行分配；
其二为后向映射法(像素填充法)，这时将输出像素逐个地映射回输入图像中，若输出像素被映射到四个输入像素之间的位置，则其灰度由它们的插值来确定。在实际中，通常采用后向映射法。
灰度级插值是用来估计像素在图像像素间某一位置处取值的过程。例如，如果用户修改了一幅图像的大小，使其包含比原始像素更多的像素，那么必须使用插值方法计算其额外像素的灰度取值。
灰度级插值的方法有很多种，但是插值操作的方式都是相同的。无论使用何种插值方法，首先都需要找到与输出图像像素相对应的输入图像点，然后再通过计算该点附近某一像素集合的权平均值来指定输出像素的灰度值。像素的权是根据像素到点的距离而定的，不同插值方法的区别就在于所考虑的像素集合不同。

最近邻插值nearest－输出像素将被指定为像素点所在位置处的像素值；
双线性插值bilinear－输出像素值是像素2×2邻域内的权平均值；
双三次插值bicubic－输出像素值是像素4×4邻域内的权平均值。

最近邻插值不够精确。双线性插值利用(x,y)点的四个最近邻像素的灰度值，按照以下方法计算(x,y)点的灰度值。

（五）步骤（计算程序流程框图）
读取原始图
建立退化图像
傅里叶变换
计算信号功率
计算噪声功率

获取噪信比k

频域分析
维纳滤波
p（m,n,）
有约束最小平方滤波方法

复原图像
复原图像
仿射变换
最近邻插值
双线型插值
三次内插法

（六）结果
（1）维纳滤波复原

图1.1 原始图像和水平运动模糊退化图像

图1.2对水平运动模糊退化图像进行维纳滤波和仿射变化

图1.3对仿射变换后的图像进行插值运算

（2）有约束最小平方滤波方法

图1.4原始图像和水平运动模糊退化图像

图1.5对约束最小平方滤波后图像进行约束最小平方滤波和仿射变化

图1.6对仿射变换后的图像进行插值运算

（七）结论与讨论
本文通过对上述三种算法 的研究和利用Matlab仿真实现得出以下结论： 从复原图像质量来看：维纳滤波复原的复原效果较差。 对于严重模糊的观测图像， 最小二乘方滤波复原可以获得较满意的复原效果，并且复原质量会更好一些。在实际应用中，要根据经验来选择最佳参数进行图像复原。对于今后图像复原算法的研究，应以提高复原算法的有效性和效率为主要研究方向，不断提高复原图像的质量和速度，并降低算法的复杂度。
维纳滤波的优点是：适用面广，无论平稳随机过程还是离散，是坐标还是向量都可以使用。对于某些难题，还可以提出滤波传递函数显式解，并进而使用由简单物理器件结合网络构成维纳滤波器。其缺点是：要得到半无限时间区间内全部观察数据前提供环境在实际中很难足够，同时它也不能应用噪声为非平稳随机过程的情况，对于响亮情况使用也不方便。
有约束最小平方滤波方法的实现过程中引入了约束条件，问题的核心是如何选择变换矩阵,选择变换矩阵的形式不同，可得到不同类型的最小平方滤波复原方法。在图像复原过程中，一般要涉及到对含噪声模糊图像的逆滤波处理，但由于系统对噪声特别敏感，使得零点附近产生震荡，并使得图像上出现过多的噪声和不必要的边缘。 因此为了减少系统对噪声特别敏感的问题，可以通过选定合适的变换矩阵以对复原后的图像强加一定程度的光滑性约束来解决，以对图像的不平滑性降至最低。

参考文献
[1]周山.MATLAB在图像复原中的应用.1.华东师范大学数学系,上海200062；2 .井冈山大学数学系，江西 吉安343009.
[2] 张亚锋,吴禄慎.图像恢复中有约束最小平方滤波技术研究.江西：南昌大学机电工程学院，2005年12月第27卷第4期
[3]于红旗，刘剑，黄知涛，周一宇.多级维纳滤波music[J].信号处理，2007年23卷6期：937-940
[4]徐飞，施晓红.MATLAB 应用图像处理[M].西安：电子科技大学出版社，2002.
[5]钟守昌1，林 洪2，赵洪洋2，彭振军2，曾仁端3.1振武汉市职工医学院 物理教研室；2.同济医科大学 附属协和医院；3.同济医科大学 物理教研室,湖北 武汉 430016.
[6]章毓晋.图像工程（上册）图像处理（第2版）.清华大学出版社.2006年.62页—76页.

附：程序代码：
clear;
close all;
%通过模拟水平运动模糊建立退化函数
d=5;
h=zeros(2*d+1,2*d+1); %分配数组
h(d+1,1:2*d+1)=1/(2*d);
% add noise
I=imread('C:\Documents and Settings\HP_USER\桌面\keshe shiyan6\mri2.jpg');
[m,n]=size(I);
fe=zeros(m+2*d,n+2*d);
fe(1:m,1:n)=I;
he=zeros(m+2*d,n+2*d);
he(1:2*d+1,1:2*d+1)=h;
F=fft2(fe);
H=fft2(he);
g=imnoise(uint8(ifft2(F.*H)),'gaussian',0,0.0001);
G=fft2(double(g));
subplot(2,4,1);imshow(I);title('原始图像');axis on;
subplot(2,4,2);imshow(g);title('水平运动模糊退化');axis on
%维纳滤波
k=0.05;
[m,n]=size(G);
for u=1:m;
for v=1:n;
i=abs(H(u,v));
i=i.^2;
s(u,v)=(1/H(u,v)*(i./(i+k)));
r(u,v)=s(u,v).*G(u,v);
end
end

r=ifft2(r);
r=uint8(real(r));
subplot(2,4,3);
imshow(r);title('维纳滤波');axis on;
%约束最小平方滤波器
[M,N]=size(G);
p=[0 -4 0;1 -4 1 ;0 1 0];
[m,n]=size(p);
pp=zeros(M,N);
pp(1:m,1:n)=p;
P=fft2(pp);
k1=0.001;
for u=1:M;
for v=1:N;
i=abs(H(u,v));
i=i.^2;
j=abs(P(u,v));
j=k1*(j.^2);
H1(u,v)=i./H(u,v);
s2(u,v)=(H1(u,v)./(i+j));
r2(u,v)=s2(u,v).*G(u,v);
end
end
r2=ifft2(r2);
r2=uint8(real(r2));
subplot(2,4,4);
imshow(r2);title('约束最小平方滤波');axis on;
%仿射变换
s=0.5;T=[s 0;0 s;0 0];
tf=maketform('affine',T);
I1=imtransform(r,tf,'bicubic','FillValues',0.3);
subplot(2,4,5);imshow(I1);title('仿射变换');axis on;

%插值
I2=imresize(I1,1,'nearest');
I3=imresize(I1,1,'bilinear');
I4=imresize(I1,1,'bicubic');
subplot(2,4,6);imshow(I2);title('最近邻插值');axis on;
subplot(2,4,7);imshow(I2);title('双线性插值');axis on;
subplot(2,4,8);imshow(I2);title('双三次插值');axis on;

tgKQcWA3PtGZ7R4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGshLs50cLmTWN60eo8Wgqv7XAv2OHUm32WGeaUwYDIAWGMeR4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGZ7R4I30kA1DkaGtgKQcWA3PtGZ7R4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGshLs50cLmTWN60eo8Wgqv7XAv2OHUm32WGeaUwYDIAWGMeR4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGZ7R4I30kA1DkaGtgKQcWA3PtGZ7R4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGshLs50cLmTWN60eo8Wgqv7XAv2OHUm32WGeaUwYDIAWGeR4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGZ7R4I30kA1DkaGtgKQcWA3PtGZ7R4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGshLs50cLmTWN60eo8Wgqv7XAv2OHUm32WGeaUwYDIAWGMeR4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGZ7R4I30kA1DkaGtgKQcWA3PtGZ7R4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGshLs50cLmTWN60eo8Wgqv7XAv2OHUm32WGeaUwYDIAWGMeR4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGZ7R4I30kA1DkaGtgKQcWA3PtGZ7R4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGshLs50cLmTWN60eo8Wgqv7XAv2OHUm32WGeaUwYDIAWGMeR4I30kA1DkaGhn3XtKknBYCUDxqA7FHYi2CHhI92tgKQcWA3PtGZ7R4I30kA1DkaG