数字图像处理-畸变校正

数字图像处理
图像畸变及校正
1 图像畸变介绍
从数字图像处理的观点来考察畸变校正, 实际上是一个图像恢复的过程, 是对一幅退化了的图像进行恢复。

在图像处理中，图像质量的改善和校正技术，也就是图像复原，当初是在处理从人造卫星发送回来的劣质图像的过程中发展、完善的。

目前，图像畸变校正的应用领域越来越广，几乎所有涉及应用扫描和成像的领域都需要畸变校正。

图像在生成和传送的过程中，很可能会产生畸变，如：偏色、模糊、几何失真、几何倾斜等等。

前几种失真主要是体现在显示器上，而后一种失真则多与图像集角度有关。

不正确的显影，打印、扫描，抓拍受反射光线的影响等方式，都会使图像产生偏色现像。

模糊、几何畸变主要是在仪器采集图片过程中产生，大多是因机器故障或操作不当影响导致，如在医学成像方面。

而几何空间失真广泛存在于各种实际工程应用中，尤其是在遥感、遥测等领域。

2 畸变产生的原因
在图像的获取或显示过程中往往会产生各种失真（畸变）：几何形状失真、灰度失真、颜色失真。

引起图像失真的原因有：成像系统的象差、畸变、带宽有限、拍摄姿态、扫描非线性、相对运动等;传感器件自身非均匀性导致响应不一致、传感器件工作状态、非均匀光照条件或点光源照明等；显示器件光电特性不一致；图像畸变的存在影响视觉效果，也是影响图像检测系统的形状检测和几何尺寸测量精度的重要因素之一。

3图像畸变校正过程所用到的重要工具
灰度直方图是关于灰度级分布的函数，是对图象中灰度级分布的统计。

灰度直方图是将数字图象中的所有像素，按照灰度值的大小，统计其所出现的频度。

通常，灰度直方图的横坐标表示灰度值，纵坐标为想像素个数。

直方图上的一个点的含义是，图像存在的等于某个灰度值的像素个数的多少。

这样通过灰度直方图就可以对图像的某些整体效果进行描述。

从数学上讲，图像的灰度直方图是图像各灰度值统计特征与图像灰度值出现的频率。

从图形上来讲，它是一个一维曲线，表征了图像的最基本的统计特征。

作为表征图像特征的信息而在图像处理中起着重要的作用。

由于直方图反映了图像的灰度分布状况，所以从对图像的观察与分析，到对图像处理结果的评价，灰度直方图都可以说是最简单、最有效的工具。

4图像颜色畸变校正介绍
图像颜色畸变现象可以是由摄像器材导致，也可以是由于真实环境本身就偏色导致，还有的是由于图像放置过久氧化、老化导致。

无论其产生的原因如何，其校正方法都是类似的。

如果用Matlab显示颜色畸变的图像RGB基色直方图，发现相对正常图像，颜色畸变的图像的直方图的三种基色的直方图中至少有一个直方图的像素明显集中集中在一处，或则集中在0处或则集中在255处，而另一部分有空缺，或则集中在中间而两边空，因此通过调整该直方图的像素点的像素值在区间[0，255]上的分布来解决图像颜色畸变问题。

如果直方图中像素集中在0一边则说明该基色偏暗，如果集中在255处则说明该基色偏亮。

下图是一有颜色畸变的图像的基色B 的直方图。

图1 基色B的直方图
很明显几乎所有像素点都集中在区间[a,b]上，这是偏暗的情况。

那么要做的是把代表基色B的矩阵的数据拉伸，使得区间[a,b]扩大为区间[a,c]。

只要做以下处理即可得到以上目的。

对每一个x ，x在[a,b]上，x*(c-a)/(b-a)，而所有的y，y 在区间[b,c]上，y=c，c=255。

MATLAB中的算法实现如下：
function dealcolor(pic,k,d)
a=imread(pic); %提取指定图像到矩阵a
b=double(a); %将矩阵a的数据转化为double型
[m,n]=size(b(:,:,k)); %取图像矩阵的行列数
fr=255/d; %设定拉伸系数
for i=1:m % 二重循环对矩阵内的每一个数据进行处理for j=1:n
bm=b(i,j,k)*fr; % 拉伸处理
if bm>255 %将所有值大于255的点都设为255
bm=255;
end
b(i,j,k)=bm;
end
end
c=uint8(b); %将矩阵b 转化为8个字节的整型数据
image(c); %显示处理过的图像
注释：a ——— 要处理的图像矩阵，是一个三维矩阵
pic ——— 要处理的图像的路径
k ——— 要处理的第几个基色矩阵，1、2、3分别代表R 、G 、B d ——— 向量，它的值是要拉伸的像素值中的最大值，0~255之间 5 图像模糊校正算法介绍
图像的灰度变化情况可以表现为一曲线。

当读入一个图像后，灰度变化就转变成了矩阵数据的变化。

反映数据变化的数学手段可以采用微分算子。

从数学的微分含义来看，“一阶微分”是描述“数据的变化率”，“二阶微分”是描述“数据变化率的变化率”。

在感应灰度变化方面，二阶微分比一阶微分更具敏感性，尤其是对斜坡渐变的细节。

因此采用二阶微分算子来处理。

最简单的各向同性微分算子是拉普拉斯微分算子。

设原图为f(x,y)，一个二维的拉普拉斯微分算子定义为：22222
y f x f f ∂∂+∂∂=∇ 将它展开就得到 )1,()1,(),1(),1(),(42+---+---=∇y x f y x f y x f y x f y x f f 。

写成图像处理运
算模版的形式就是⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡----=010141010L 。

设处理后的图像为g(x,y)，则 ),(),(),(2y x f y x f y x g ∇+=。

用模板表示则是⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡----=010*******L ，模板中心点就是要处理的像素点。

如果连对角线方向都考虑的话，模板的形式就表现为⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡--------=1111911111L 。

由于模板是一个3阶矩阵，所以模板处理不了图像矩阵边缘的点。

将这些处理不了的点像素值都设置为255。

在Matlab 中的算法实现为：
function pic_tem = pick_tem(b,x,y,j,tem_n) % b 是一个三维矩阵，j 是代表第几层矩阵
%tem_n代表模板号0代表采用第一种模板，1代表采用第二种模板
if tem_n==0 %第一种模板
%返回二阶微分处理后的值
pic_tem=5*b(y,x,j)-(b(y,x-1,j)+b(y+1,x,j)+b(y-1,x,j)+b(y,x+1,j));
else %第二种模板
%返回二阶微分处理后的值
pic_tem=9*b(y,x,j)-(b(y-1,x-1,j)+b(y-1,x,j)+b(y-1,x+1,j)+b(y,x-1,j)+b(y+1,x+1,j)+b(y +1,x-1,j)+b(y+1,x,j)+b(y+1,x+1,j));
end
function faintness(pic,tem_n)
%pic代表处理的图片的路径，tem_n代表采用第几种模板
a=imread(pic); %导入图像，并把数据存放到三维矩阵a中
b=double(a); %将数据转化为双精度
[h,w,k]=size(b); %取矩阵的维数
g=zeros(h,w,k)+255; %生成一个与原图像矩阵有相同维数的矩阵g
for y=2:h-1 %剔除图像矩阵最外一层的点，处理可以被模板包含的点
for x=2:w-1
for j=1:k
g(y,x,j)=pick_tem(b,x,y,j,tem_n); %采用二阶微分算子处理，返回处理点（x,y）后的像素值end
end
end
image(uint8(g));%显示处理后的图像
6图像几何畸变校正算法介绍
对于畸变的光学系统，畸变空间中的直线在像空间中一般不再是直线，而只有通过对称中心的直线是例外。

因此在进行桶形畸变校正时须先找出对称中心，再进行通用的几何畸变校正过程。

桶形畸变校正一般步骤：
（1）找出畸变图对称中心，将畸变图代表的地址空间关系转换为以对称中心为原点的空间关系。

（2）空间变换：对输入图像（畸变图）上像素重新排列以恢复原空间关系。

即利用地址映射关系为校正图空间上的每一个点找到它们在畸变图空间上的对应点。

（3）灰度插值：对空间变换后的像素赋予相应的灰度值以恢复原位置的灰度值。

几何畸变的校正要使用几何（坐标）变换，包括平行移动、旋转、扩大缩小等简单的变换。

在这里，先在直角坐标系下研究地址映射关系，在程序中则使用矩阵坐标系。

设[f(x,y)]是原图，[f(u,v)]是发生畸变后的图像。

畸变校正的基本思想是，找出由(u,v) (x,y)的坐标变换αT (地址映射),然后令：)),((),(v u T f y x f α= α为参数向量。

若畸变只是简单的纵横比的改变和倾斜，那么仿射变换可以校正这种畸变。

取：{}F E D C B A T ,,,,,:=αα
[]C Bx Au x ++= ][F Ev Du y ++= []x 表示去最接近x 的整数
若能得到参数α的估计，问题就可以解决了。

但一般的畸变都不只是简单的纵横比的改变和倾斜，通常遇到的都是空间扭曲型几何畸变，俗语讲就是橡胶层面拉伸。

它是曲线畸变，这里采用二次多项式来模拟它,可表示为：
25423210y a xy a x a y a x a a u +++++=(1) 25423210y b xy b x b y b x b b v +++++=(2) 同样，只要可以取得参数),(i i b a 的估计，畸变函数便可知，那么原则上，可以通过上式多项式变换来获得修正的空间扭曲映射。

(1)和(2)都是一个有6个参数的二元二次方程，那么只要在畸变图和校正图上各取六对对应点（其中从校正图上的点是估计值）就可以通过解方程组得到参数),(i i b a 的估计，理论上，取的对应点对数越多得到的参数),(i i b a 估计就越精确。

设取m 对对应点，用向量来表示为
],...,,[21m t u u u U = ],...,,[21m t v v v V =
⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=2222222222
211121111.....
.......11m m m m m m y y x x y x y y x x y x y y x x y x A 其系数 ],...,,[10m t a a a a = ],...,,[10m t b b b b =
假设所取的m 对对应点组成的矩阵A 可逆，也就是说m 对对应点线性不相关，则容易计算得到系数),(i i b a ：U A a -= V A b -=
由于采用灰度插值，所以在校正的处理过程是对校正图上每一点映射到畸变图，然后通过灰度插值来得到这一点的灰度值。

因此校正所采用的是地址逆映射：
[]a y xy x y x u *,,,,,122= []
b y xy x y x v *,,,,,122= (∆)
由地址映射(∆)计算得到的（u,v ）可能是非整数，而畸变图[f(u,v)]是数字图像，其像素值仅在坐标为整数处有定义，所以在非整数处的像素值要用其周围一些整数处的像素值来计算，这叫灰度插值。

灰度插值有邻近插值法、双线性插值法等等。

邻近插值法得到的图像通常都会出现锯齿现象，而双线性插值法得到的图像比邻近插值法得到的图像精确得多。

通常情况下，双线性插值的精确度已经可以满足一般图像处理的要求，并不需要更高精确度的灰度插值。

因此在这里将采用双线性插值法。

双线性插值利用(u,v)周围的四个最邻近像素的灰度值，根据下面方法来计算(u,v)处的灰度值。

设(u,v)四个邻近像素点为ABCD 坐标分别为(i,j), (i+1,j), (i,j+1),(i+1,j+1)
设 i u -=α j v -=β
首先先计算出E,F 两处的灰度值，f(E)和f(F)
f(E)=β[f(C)-f(A)]+f(A) f(F)=β[f(D)-f(B)]+f(B)
再计算(u,v)： f(u,v)=α[f(F)-f(E)]+f(E)
此f(u,v)值代表的就是校正后图像中(x,y)处的灰度值。

算法的流程为：
B C
D E F
(i,j) (i,j+1) (i+1,j)
在Matlab中的算法实现为：
function gmodify(pic,uv,gm,og) %pic表示要处理的图像的路径文件名
%uv是一个二维矩阵,uv(:,1)代表上面提到的t U,uv(:,2)表示t V
%gm是一个二维矩阵,gm(j,:)代表在校正图空间上与uv(j,:)一一应的点
%og 代表对称中心，它是一个二维向量
a=imread(pic);
b=double(a);
n=size(gm(:,1));
for k=1:n%转换到以对称点为原点的空间关系并构造矩阵A
A(k，：)=[1,gm(k,1)-og(1),gm(k,2)-og(2),gm(k,1)-og(1)^2, (gm(k,1)-og(1))*(gm(k,2)-og(2)),(gm(k,2)-og(2))^2];
end
[h,w]=size(b(:,:,1));
sp=zeros(h,w,3)+255;
a0=pinv(A)* uv(:,2); %计算上面提到的地址映射的系数估计a
b0=pinv(A)* uv(:,1); %计算上面中提到的地址映射的系数估计b
for i=1:h %从理想图像矩阵出发处理
for j=1:w
x=[1,j-og(1),i-og(2),(j-og(1))^2,(i-og(2))*(j-og(1)),(i-og(2))^2];
u=x*a0+og(2); % 逆向映射（j，i）到畸变图像矩阵（v，u）
v=x*b0+og(1);
if (u>1)&&(u<w)&&(v>1)&&(v<h) %处理在图像大小范围内的像素点
uu=floor(u); %对u取整
vv=floor(v); %对v取整
arf=u-uu; %计算上面提到的α
bta=v-vv; %计算上面提到的β
for k=1:3 %进行灰度双线性插值
ft1=(1-bta)*b(vv,uu,k)+bta*b(vv+1,uu,k);
ft2=(1-bta)*b(vv,uu+1,k)+bta*b(vv+1,uu+1,k);
sp(i,j,k)=(1-arf)*ft1+arf*ft2;
end
end
end
end
image(uint8(sp)); %显示校正图像
7图像几何畸变校正实例
下面是一幅已经发生严重几何变形的图片，根据上面描述的算法，对该图片进行几何畸变校正。

图8 畸变图
这里选择网状图是因为网状图更容易分析，校正效果也更容易观察。

观察畸变图，发现从上数下第五条横线和从左数起第六条竖线都是直线，因此选用这两条直线的交叉点作为对称中心进行畸变校正。

然后，从畸变图中选出4个线性不相关的点，并估计它们在校正图上的对应点。

因此得到的三个输入参数的值分别如下：
pic=’tt.jpg’;
uv=[144,26;26,198;144,416;308,241;19,29;303,30;307,325];
gm=[144,29;29,198;144,412;304,241;15,26;306,26;307,327];
og=[145，241];
gmodify(pic,uv,gm,og);
得到的校正图如下：
图 9 校正后的图片
虽然经过校正得到的图像仍有一定程度的畸变，但是畸变程度小多了。

对几何畸变图像的校正计算, 不可能完全准确的恢复原来的图像。

这是因为影响畸变校正精度的因素是多方面的。

其中对校正算法精度影响较大的几种误差有: 对称中心估算误差、系数估计误差、插值误差和灰度校正误差等。

该校正方法关键在于对称中心的估计和地址映射的系数估计，而系数估计的关键是在于所取对应点对的精确程度。

显然，校正后的图像仍有一定程度的畸变是由于对应点对不够精确。

校正图上的对应点的估计是很困难的，通常都需要对程序运行的结果图进行分析，然后不断进行调整。

。