摘要:光学成像系统非线性几何畸变的高精度数字校正仍然是一个未能很好解决的问题。其中,衡量畸变程度的参数难以精确测量是最重要的原因之一。在以径向几何畸变为主的非线性几何畸变模型中,通过对影响畸变参数测量精度的各种因素的分析,提出了提高畸变参数测量精度的方法。详细介绍了通过计算机自动测量畸变参数的算法,并给出了实现数字校正的算法。实验表明,能够比较精确地测出实现畸变校正所需的各参数。应用到不规则平面物体面积的测量中,获得了很好的效果。
1 引 言
图像畸变广泛存在于应用广角镜头的成像系统中。畸变的存在不利于图像的辨认、分析和判断。在一些定量分析领域,需要对畸变进行校正,则校正精度将直接影响定量分析的精度。例如,需要测量内窥镜图像上病变区域的大小时,需要依靠机器视觉进行精确定位时[1],测量不规则平面物体的面积时[2],等等。成像系统的几何畸变分为两类,线性几何畸变和非线性几何畸变。线性几何畸变的校正已经得到了比较好的解决,但是非线性几何畸变的校正还是一个未能成功解决的课题[3]。
为了解决这一问题,国内外很多学者进行了这方面的研究[1,3～6]。从这些文献可以看出,成像系统对称中心(以下简称对称中心)的确定是实现畸变校正的前提。在确定对称中心的方法中,有的采用点阵图案的校正样板[1,4],根据样板上黑点的中心找到三条最接近于直线排列的点列,用三条直线拟合相应的点列,以这三条直线的交点作为对称中心[1]。由于各点列之间必然存在一定的距离,所确定的对称中心必存在较大误差。有的采用栅格状的校正样板,根据样板上经纬线的弯曲程度估计对称中心位置[5]。同样由于各经纬线之间存在距离,难以得到精确的对称中心。有的利用平行栅[3,6],根据水平放置时图像上某一条曲线的顶点确定对称中心的水平位置,类似的,将平行栅垂直放置,再确定其垂直位置[3]。由于曲线在顶点处变化缓慢,必存在较大误差。而且,当对称中心与光敏面的中心不重合时,单凭主观观察很难保证平行栅水平或垂直,需借助于专门的仪器。
我们提出一种新的方案来提高对称中心的测量精度,进而提高畸变校正的精度。利用倾斜放置的平行栅,在任意的两个不同的倾斜角下,得到两幅存在畸变的图像。对每一幅图像,找到最接近于过对称中心的曲线,以这两条曲线的交点作为对称中心。然后,在以径向几何畸变为主的非线性几何畸变模型中,用最小二乘法求出系统畸变参数。当对称中心确定后,再精确测定图像采集控制器的纵、横采集比。这种方法无须借助于专门的仪器,不存在主观误差,更重要的是,可以比较准确地测出实现畸变校正所需的各个参数。将所得参数应用到不规则平面物体面积的测量中进行畸变校正,测量精度得到了很大提高。
2 几何畸变数字校正原理
2.1 成像系统简介
如图1所示,成像系统由摄像头、图像采集控制器、计算机、显示器等几部分组成。摄像头将成像到其光敏面上的物体的像转换成电信号,经由图像采集控制器,显示到显示器上。图像采集控制器可以将物体的像采集下来,以数字形式保存到存贮介质,也可以将存贮介质中数字化的图像在显示器上显示出来。
2.2 镜头畸变
镜头的几何畸变有径向和切向两类,通常径向畸变比切向畸变大得多。在以径向几何畸变为主的非线性几何畸变模型中,无畸变像面坐标和畸变像面坐标之间的关系为[3]
其中(x,y)为无畸变像平面坐标;(xd,yd)为畸变像平面坐标,其坐标原点都在对称中心;k为径向畸变系数。
2.3 图像采集控制器畸变
图像采集控制器畸变主要是纵、横方向的采集间距不相等引起的。引入图像平面坐标系(xf,yf),以及纵、横方向坐标变换系数sy、sx。假设摄像头的光敏面和像平面重合。但是,光敏面的中心与对称中心一般并不重合。设对称中心在图像平面坐标系中的坐标为(cx,cy)。在图像平面坐标系中,存在镜头畸变的坐标(xd,yd)将进一步产生畸变,成为(x′f,y′f),满足[3]
2.4 有畸变和无畸变图像坐标关系
如果没有图像采集控制器畸变,即sx=sy,则无畸变像平面上一点(x,y)在图像采集器平面坐标系中的坐标(xf,yf)应该为
将(2)式、(3)式代入(1)式可得有畸变和无畸变图像之间的坐标关系为(文献[3]中的(4)式和(6)式存在错误)
其中Δ=cy(1 -sy/sx);K=k/s2x;K称为系统畸变参数。为了方便,以后将(cx,cy)、K以及sy/sx统称为畸变参数。
2.5 畸变校正
当已知各个畸变参数时,即可进行畸变校正。畸变校正可以分为两步,第一步是对理想图像所在的坐标空间按照(4)式进行几何变换,使理想图像上的像素点与实际图像上的点对应起来,即给定(xf,yf),求出对应的(x′f,y′f).第二步是进行灰度校正。因为变换后的点不可能总是落在实际图像的像素点上,所以需要根据实际图像上该点附近各像素点的灰度值来估计该点的灰度值,从而得到与之对应的理想图像上像素点的灰度值,称为灰度校正。我们采用二元全区间插值法[7]进行灰度校正,方法如下:
设实际图像上各点的灰度用z(x′,y′)表示,以变换后的点(x′f,y′f)为中心,在X方向和Y方向前后各取四个坐标
3 畸变参数的测量步骤和算法
3.1 测量(cx,cy)的步骤和算法
首先让平行栅在两个不同倾斜角下拍摄两幅图像,然后对每一幅图像按照如下步骤求取一条曲线:
1)用平均灰度作为阈值进行二值化。
2)以不同的编号,将各条曲线区分开来。
3)垂直方向压缩各曲线。所用平行栅线条宽度为几个像素,可以保证不至于将同一条曲线分割开来。垂直方向压缩后,可以减小以后的计算量,也可以减小因取整造成的纵坐标误差。
4)判断各曲线的弯曲程度。对每一条曲线,可以按照最小二乘法拟合一条直线,设可以表示为y=aix+bi,则可以用
表示该曲线的弯曲程度。
5)根据第4)步结果,找到最小的δi,取其对应曲线两侧的曲线,设如图2中的实线所示,它们的弯曲程度分别为δ1、δ2.
6)利用如下公式产生第三条曲线上各点的坐标
如图2中虚线所示。计算其弯曲程度δ3。若δ1<δ2,则用第三条曲线取代第二条曲线,反之,则取代第一条直线。
7)重复第6)步,直到可以认为两条曲线已经重合。则该重合曲线非常接近于通过对称中心。由于两幅图像倾斜角不同,经过以上步骤所确定的两条曲线必然相交,求出其交点,即可作为对称中心。
3.2 测量K的步骤和算法
理论上,知道了(cx,cy)、sy/sx之后,只要在存在畸变的平行栅图像上的同一条曲线上取三点(x′fi,y′fi),i= 1,2,3,由(4)式解出(xfi,yfi),并让它们位于同一条直线上,即可解出K。但是,即便能够给出该复杂方程,通过这么少的点求出的K也必然存在很大误差。
我们采用最小二乘法确定K。方法是,先假定一个K,然后对存在畸变的平行栅图像进行校正,求取各曲线的总弯曲程度。尝试改变K,直到各曲线的总弯曲程度减小到最小。
这样做似乎需要先知道sy/sx。但是,由于图像采集控制器畸变是线性的,像平面上的直线经镜头畸变后仍然是直线,所以,在采用上述方法求K时,完全可以先假定sy/sx= 1。此时,(4)式变为
这样做,不仅可以减小计算量,还可以减小计算误差对K的影响。具体步骤和算法如下:
让平行栅近于水平放置,以便得到尽量长的曲线,拍摄一幅图像。对其进行如下处理：
1)假设一个K。按照(8)式和(5)式对图像进行校正。
2)对校正后的图像进行同3.1中的1)～4)步的处理。
3)按照下式计算各曲线总弯曲程度
4)采用变步长法搜索K,求出使δ最小的K作为系统畸变参数。
3.3 测量sy/sx的方法
当已知(cx,cy)之后,可以用一组同心圆,使圆心位于(cx,cy)附近,拍摄一幅图像。对所得图像进行搜索,找到同一个圆上最上边点的纵坐标和最下边点的纵坐标,用二者的差值除以该圆上最右边点的横坐标与最左边点的横坐标的差,作为由该圆确定的sy/sx。对所有圆求sy/sx的平均值作为图像采集控制器的纵、横采集比。
4 实验结果
依据前述算法,我们首先用Visual C++ 5.0作为开发工具,编写了基于Windows 95/98的测量畸变参数和进行畸变校正的应用程序。然后进行了实际测量。
实验装置为:摄像头型号为WV-CP410/G,图像采集控制器为CA6300图像卡,分辨率选为768×576,计算机为奔腾133,镜头为f=3.5 mm、67°的短焦距广角镜头。照明条件为室内日光灯加自然光。平行栅和同心圆都是自己制备的,用CorelDRAW 6.0绘制,用激光打印机打印,分别如图3(a)、(b)所示。平行栅线条间距约1.32 mm,同心圆半径以约1.25 mm递增。实验前,调整物平面尽量与成像系统光轴垂直,实验过程中像距物距保持不变。
首先让平行栅从左上至右下倾斜,在不同的角度下拍摄了4幅图像,分别命名为a1～a4。然后让其从左下至右上倾斜,在不同的角度下拍摄了4幅图像,分别命名为b1～b4。又让其在水平位置附近,拍摄了8幅图像,分别命名为c1～c8。在前两组中任意组合选出两幅图像,以图像左上角为原点,计算(cx,cy),在此基础上,从第三组中顺序选一幅图像,计算K,结果如表1所示。
知道了(cx,cy)之后,让同心圆圆心在其附近,拍摄了10幅图像,分别命名为d1～d10。对每幅图像计算sy/sx,结果如表2所示。然后,又取各畸变参数的测量平均值对这10幅图像进行了校正,并对校正后的图像重测了纵横采集比,结果在表2中用(sy/sx)′同时给出,可以看出,已经非常接近于1。另外,还让同心圆圆心在图像中心位置附近拍摄了一幅图像,圆心位置为(383,282),用它测得sy/sx= 1.013 6,和平均值的差值为0.001 4,几乎是标准方差σs的3倍。但是,按照各畸变参数的测量平均值对其进行了校正后再测(sy/sx)′,结果为0.999 0。
我们认为,可以通过计算被测量的标准方差与其平均值的相对大小来估计测量结果的好坏。对K来说,δK=σK/K= 1.0%。对cx、cy来说,应将坐标原点平移到图像中心计算,结果分别为:
同样,对sy/sx来说,应为δs=σs/ 1 - (sy/sx) =3.3%。究其原因,可能是求K时参加运算的点最多,使得δK最小,求sy/sx时参加运算的点最少,使得δs最大。求cx、cy时对曲线进行了垂直方向的压缩,减小了纵坐标因取整效应带来的误差,使得δcy小于δcx。考虑到可能存在其它各种线性、非线性畸变,以及各种计算误差,可以说测量结果具有较高的精度。
为了从直观上观察畸变校正的效果,又制备了栅格样板,其畸变图像如图4(a)所示,校正后的图像如图4(b)所示。可见校正效果很好。只是图像边缘有少量信息丢失,因为矫正后的图像大小仍然是768×576。如果采用比较大的分辨率,则会在边缘出现部分像素无图像信息。这是畸变矫正不可避免的。
为了从数值上说明矫正效果,从图4(b)中选取了11个完整的方格(右下角的不完整,未选),它们是第1、3、5行的第1、3、4、6个方格,以及图4(a)中对应的11个方格,按照行列号编号为Rij。在表3中给出了它们的45°、135°方向对角线相对于图4(b)中方格R34的45°对角线的相对误差,以及变形最大的纵、横边相对于图4(b)中方格R34的纵边的相对误差。各线长度是在PhotoShop4.0下按照400%显示后,用其画线工具在两顶点之间画直线而得到。矫正后误差大大减小,由原来的最大18%降低为最大0.9%。
将所测得的畸变参数应用到不规则平面物体面积的测量中,取得了很好的效果。测量原理参见文献[2]。用一个长41.3 mm,宽30.1 mm的长方形作为被测物体,先让其在图像中心位置附近拍摄一幅图像,命名为e0,然后在其它位置,包括旋转不同的角度,拍摄了10幅图像,分别命名为e1～e10。用e0作为标准对系统进行标定,对存在畸变的e1～e10的测量结果如表4所示。可见测量误差比较大。然后利用测得的畸变参数对所有图像进行畸变校正,用校正后的e0作为标准对系统进行标定,测量校正后的e1～e10,结果也在表4中给出。可见,测量精度得到了很大提高。
5 结 论
利用我们提出的新方案来确定成像系统对称中心的位置,可以提高测量精度,进而可以提高其它畸变参数的测量精度,从而可以提高畸变校正的精度。而且,这种方法对校正样板的放置位置要求不严格,即使不借助于专门的仪器,也不存在主观误差。应用到不规则平面物体面积的测量中,测量精度得到了很大提高。
参考文献:
[1]刘航,郁道银,杜吉等.广角成像系统光学畸变的数字校正方法[J].光学学报,1998,18(8):1108～1112.
[2]张全法,韩要轩,杨海彬等.用面阵CCD测量不规则平面物体的面积[J].仪表技术与传感器,2000,(2):36～39.
[3]林子扬,王继清,黄乐天.一种光学成像垂轴几何畸变校正方法[J].光学技术,1998,(5):19～22.
[4] Smith W E, Vakil N, Maislin S A. Correction of distortion in en-doscope images [J]. IEEE Trans. Med. Imag., 1992, MI-11(1):117～122.
[5] Hideaki Haneishi, Yutaka Yagihashi, Yoichi Miyake. A newmethod for distortion correction of electroNIc endoscope images[J]. IEEE Trans. Med. Imag. , 1995, MI-14(3):548～555.
[6] Yoshihiko Nomura, Michihiro Sagara, Hiroshi Naruse,et al. Sim-ple Calibration algorithm for high-distortion-lens camera [J].IEEE Trans. Pattern Anal. &Matine Intell. , 1992, PAMI-14(11):1095～1099.
[7]徐士良. C常用算法程序集(第二版)[M].北京:清华大学出版社,1996.
作者简介:张全法(1966-),男,河南人,郑州大学物理工程学院讲师,硕士,从事传感器应用技术、计算机应用技术研究。
