第3章像素空间关系 3.1像素间的基本关系 图像的基本组成单元是像素,像素在图像空间中按照某种规律排列,有一定的相互联系。常见的像素间的基本关系包括像素的邻域、邻接、连接、通路和连通以及像素集合间的邻接、连接和连通。 3.1.1像素的邻域 像素的邻域是指一个像素的相邻像素构成的像素集。邻域的类型一般包括4邻域、对角邻域和8邻域,如图31所示。 图31邻域的类型 1. 4邻域 像素p的4邻域是指由像素p的水平方向(即左右)和垂直方向(即上下)共4个相邻像素所组成的集合,记为N4(p)={r1,r2,r3,r4}。若像素p的坐标为(x,y),则像素p的4邻域像素的坐标分别为r1: (x-1,y)、r2: (x,y-1)、r3: (x+1,y)和r4: (x,y+1),如图31(a)所示。 2. 对角邻域 像素p的对角邻域是指由像素p的对角方向(即左上、左下、右上、右下)共4个相邻像素所组成的集合,记为ND(p)={s1,s2,s3,s4}。若像素p的坐标为(x,y),则像素p的对角邻域像素的坐标分别为s1: (x-1,y+1)、s2: (x-1,y-1)、s3: (x+1,y-1)和s4: (x+1,y+1),如图31(b)所示。 3. 8邻域 像素p的8邻域是指由像素p的水平、垂直和对角方向(即左右、上下、左上、左下、右上、右下)共8个相邻像素所组成的集合,记为N8(p)={r1,r2,r3,r4,s1,s2,s3,s4}。若像素p的坐标为(x,y),则像素p的8邻域像素的坐标分别为r1: (x-1,y)、r2: (x,y-1)、r3: (x+1,y)、r4: (x,y+1)、s1: (x-1,y+1)、s2: (x-1,y-1)、s3: (x+1,y-1)和s4: (x+1,y+1),如图31(c)所示。 需要注意的是,如果像素p本身处于图像的边缘,则它的4邻域N4(p),对角邻域ND(p)和8邻域N8(p)中的若干个像素将位于图像之外。 3.1.2像素的邻接 像素的邻接是指一个像素与其邻域中的像素的接触关系。邻接的类型根据邻域的类型的不同一般分为4邻接、对角邻接和8邻接。 1. 4邻接 4邻接是指一个像素与其4邻域中的像素的接触关系。若两个像素为p和r,则像素p与像素r满足4邻接可表示为 p∈N4(r)(31) 注意: 像素p与像素r满足4邻接等价于像素r与像素p满足4邻接,即式(32)成立。 p∈N4(r)r∈N4(p)(32) 2. 对角邻接 对角邻接是指一个像素与其对角邻域中的像素的接触关系。若两个像素为p和r,则像素p与像素r满足对角邻接可表示为 p∈ND(r)(33) 注意: 像素p与像素r满足对角邻接等价于像素r与像素p满足对角邻接,即式(34)成立。 p∈ND(r)r∈ND(p)(34) 3. 8邻接 8邻接是指一个像素与其8邻域中的像素的接触关系。若两个像素为p和r,则像素p与像素r满足8邻接可表示为 p∈N8(r)(35) 注意: 像素p与像素r满足8邻接等价于像素r与像素p满足8邻接,即式(36)成立。 p∈N8(r)r∈N8(p)(36) 需要注意的是,邻接仅考虑了像素间的空间关系,与像素的属性值无关。 3.1.3像素的连接 两个像素的连接是指两个像素必须邻接(即接触)且它们的属性值必须满足某个特定的相似准则。属性值一般采用像素的灰度值。相似准则可以是灰度值相等,或者同在一个灰度值集合中取值,记为V。例如,在一张二值图像中,定义两个灰度值为1的像素之间的连接,可以取相似准则为灰度值集合V={1}; 在一张256色的灰度图像中,定义灰度值为100~105的像素之间的连接,可以取相似准则为灰度值集合V={100,101,102,103,104,105}。 连接的类型根据邻域的类型的不同一般分为4连接、对角连接、8连接以及混合连接。 1. 4连接 4连接是指两个像素4邻接且它们的属性值满足某个特定的相似准则。若两个像素为p和r,像素的属性值函数为f(),相似准则为V,则4连接的条件可表示为 (p∈N4(r))∧(f(p)∈V)∧(f(r)∈V)(37) 2. 对角连接 对角连接是指两个像素对角邻接且它们的属性值满足某个特定的相似准则。若两个像素为p和r,像素的属性值函数为f(),相似准则为V,则对角连接的条件可表示为 (p∈ND(r))∧(f(p)∈V)∧(f(r)∈V)(38) 3. 8连接 8连接是指两个像素8邻接且它们的属性值满足某个特定的相似准则。若两个像素为p和r,像素的属性值函数为f(),相似准则为V,则8连接的条件可表示为 (p∈N8(r))∧(f(p)∈V)∧(f(r)∈V)(39) 4. 混合连接 混合连接又称m连接,是指两个像素的属性值必须满足某个特定的相似准则且满足下列两个条件之一: ①两个像素4邻接; ②两个像素对角邻接且它们4邻域的交集在相似准则的意义下是空集。若两个像素为p和r,像素的属性值函数为f(),相似准则为V,则混合连接的条件可表示为式(310)或式(311): (p∈N4(r))∧(f(p)∈V)∧(f(r)∈V)(310) (p∈ND(r))∧(f(ND(p)∩ND(r))∩V=)(311) 例3.1混合连接的判定。一张4×4的二值图像模板如图32(a)所示,相似准则V={1},当图像数据如图32(b)和图32(c)取值时,判定像素p和r是否满足混合连接的条件。 图32混合连接的判定 解: 由图像模板可知,像素p和r互为对角邻接,即p∈ND(r),且像素p和r对角邻域的交集为: ND(p)∩ND(r)={c,d}。 (1) 当图像数据如图32(b)取值时,有: f(c)=0,f(d)=0,则 f(ND(p)∩ND(r))=f({c,d})=f(c)∪f(d)={0}∪{0}={0} 因此,f(ND(p)∩ND(r))∩V={0}∩{1}= 所以,像素p和r满足混合连接的条件。 (2) 当图像数据如图32(c)取值时,有: f(c)=1,f(d)=0,则 f(ND(p)∩ND(r))=f({c,d})=f(c)∪f(d)={1}∪{0}={0,1} 因此,f(ND(p)∩ND(r))∩V={0,1}∩{1}={1}≠ 所以,像素p和r不满足混合连接的条件。 混合连接可以认为是8连接的一种变型,引进混合连接是为了消除使用8连接时常出现的多路问题。 例3.2多路问题示例。一张3×3的二值图像模板如图33(a)所示,相似准则V={1},像素a、b、c、d的灰度值均取1,像素e的灰度值取0。试分析采用8连接和混合连接时像素a、b、c、d的连接情况。 图33多路问题示例 解: 采用8连接时,像素a、b、c、d的连接情况如图33(b)所示。此时,像素a和c之间存在两条8连接通路,即abc和ac,导致像素a和c之间产生多路问题。 采用混合连接时,像素a、b、c、d的连接情况如图33(c)所示。此时,像素a和b之间满足混合连接条件,像素b和c之间满足混合连接条件,因此,像素a和c之间存在一条混合连接通路abc。 另外,f(ND(a)∩ND(c))=f({b,e})=f(b)∪f(e)={1}∪{0}={0,1} 因此,f(ND(a)∩ND(c))∩V={0,1}∩{1}={1}≠ 所以,像素a和c之间不满足混合连接条件,即不存在混合连接通路ac。 因此,采用混合连接时,像素a和c之间不存在多路问题。 3.1.4像素的通路 从一个具有坐标(x,y)的像素p到另一个具有坐标(s,t)的像素q的一条通路定义为由一系列具有坐标(x0,y0),(x1,y1),…,(xn,yn)的独立像素组成的集合,并且满足以下条件: (1) (x0,y0)=(x,y),(xn,yn)=(s,t); (2) (xi,yi)与(xi-1,yi-1)邻接; (3) 1≤i≤n,n为通路的长度。 根据邻接类型的不同,通路的类型也不同。如果像素间邻接的类型为4邻接、对角邻接或8邻接,则对应通路的类型为4通路、对角通路或8通路。 3.1.5像素的连通 像素的连通是指像素通路上所有像素的属性值均满足某个特定的相似准则。根据连接类型的不同,连通的类型也不同。如果像素间连接的类型为4连接、对角连接、8连接或混合连接,则对应连通的类型为4连通、对角连通、8连通或混合连通。 像素连接可以看作是像素连通的一种特例。当n=1时,两个连通的像素也是连接的。 3.1.6像素集合的邻接 图像可以看作是像素的集合。像素集合根据一定的规则可以划分成若干个子集,每一个子集就构成了一个子图像。对于两个图像子集S和T来说,如果S中的一个或一些像素与T中的一个或一些像素邻接,则称两个图像子集S和T是邻接的。如果像素间邻接的类型为4邻接、对角邻接或8邻接,则对应图像子集的邻接类型为4邻接、对角邻接或8邻接。 3.1.7像素集合的连接 两个图像子集S和T是连接的是指两个图像子集必须邻接且邻接像素的属性值必须满足某个特定的相似准则。也就是说,对于两个图像子集S和T来说,如果S中的一个或一些像素与T中的一个或一些像素连接,则称两个图像子集S和T是连接的。如果像素间连接的类型为4连接、对角连接、8连接或混合连接,则对应图像子集的连接类型为4连接、对角连接、8连接或混合连接。 3.1.8像素集合的连通 对于图像子集S中的两个像素p和q来说,如果存在一条完全由S中的像素组成的从p到q的通路,则称p在S中与q连通。对于S中的任一个像素p,所有与p相连通且又在S中的像素的集合(包括p)称为S中的一个连通组元。图像里同一个连通组元中的任意两个像素相互连通,而不同连通组元中的像素互不连通。 图像中每个连通组元构成图像中的一个区域。图像可以认为是由一系列区域组成。区域的边界也称区域的轮廓,是该区域的一个子集,它将该区域与其他区域分开。组成一个区域边界的像素本身属于该区域而在其邻域中存在不属于该区域的像素。 3.2像素间的距离 像素间的距离是指像素在空间的接近程度。设3个像素为p、q和r,坐标分别为(x,y)、(s,t)和(u,v),则距离量度函数D必须满足下列三个条件: (1) D(p,q)≥0; (2) D(p,q)=D(q,p); (3) D(p,r)≤D(p,q)+D(q,r)。 其中,条件(1)表明两个像素之间的距离总是为正值,若D(p,q)=0当且仅当p=q; 条件(2)表明像素之间的距离与起点和终点的选择无关; 条件(3)表明像素之间的最短距离是沿直线的。 常见的距离度量方法包括欧氏距离、城区距离、棋盘距离和混合距离。 3.2.1欧氏距离 欧氏(Euclidean)距离记为DE,设像素点p的坐标为(x,y),像素点q的坐标为(s,t),则像素点p和q之间的欧氏距离的定义如下式所示: DE(p,q)=(x-s)2+(y-t)2(312) 欧氏距离的几何意义为: 距离坐标为(x,y)的像素的DE距离小于或等于某个值d的像素都包括在以(x,y)为中心以d为半径的圆中,如图34所示。 图34欧氏距离 其中,图34(a)为距离坐标为(x,y)的像素的DE距离小于或等于3的像素所构成的圆形区域,其值为各像素点距离中心像素点的距离值,该值已经过四舍五入处理; 图34(b)为距离值对应的3D透视图。 3.2.2城区距离 城区(cityblock)距离记为D4,设像素点p的坐标为(x,y),像素点q的坐标为(s,t),则像素点p和q之间的城区距离的定义如下式所示: D4(p,q)=|x-s|+|y-t|(313) 城区距离的几何意义为: 距离坐标为(x,y)的像素的D4距离小于或等于某个值d的像素都包括在以(x,y)为中心的菱形中,如图35所示。其中,图35(a)为距离坐标为(x,y)的像素的D4距离小于或等于3的像素所构成的菱形区域,其值为各像素点距离中心像素点的距离值; 图35(b)为距离值对应的3D透视图。 图35城区距离 由城区距离的定义可知,距离像素p的城区距离为1的像素就是像素p的4邻域像素。因此,像素p的4邻域可以通过城区距离定义为 N4(p)={r|D4(p,r)=1}(314) 式中,r为某个像素。 3.2.3棋盘距离 棋盘(chessboard)距离记为D8,设像素点p的坐标为(x,y),像素点q的坐标为(s,t),则像素点p和q之间的棋盘距离的定义如下式所示: D8(p,q)=max(|x-s|,|y-t|)(315) 棋盘距离的几何意义为: 距离坐标为(x,y)的像素的D8距离小于或等于某个值d的像素都包括在以(x,y)为中心的正方形中,如图36所示。 图36棋盘距离 其中,图36(a)为距离坐标为(x,y)的像素的D8距离小于或等于3的像素所构成的正方形区域,其值为各像素点距离中心像素点的距离值; 图36(b)为距离值对应的3D透视图。 由棋盘距离的定义可知,距离像素p的棋盘距离为1的像素就是像素p的8邻域像素。因此,像素p的8邻域可以通过棋盘距离定义为式(316)。 N8(p)={r|D8(p,r)=1}(316) 式中,r为某个像素。 3.2.4混合距离 像素点之间的欧氏距离(DE)、城区距离(D4)和棋盘距离(D8)刻画了像素在空间的接近程度,其大小仅与像素的坐标有关,与像素本身及其邻近像素的属性值无关。两个像素点p和q之间的欧氏距离DE(p,q)等于它们之间的直线距离,通常无法对应某条具体通路。城区距离D4(p,q)等于它们之间最短的4通路的长度。棋盘距离D8(p,q)等于它们之间最短的8通路的长度。像素间距离与通路示意图如图37所示。 图37像素间距离与通路示意图 像素点之间的混合距离Dm同样刻画了像素在空间的接近程度,其大小不仅与像素的坐标有关,还与像素本身及其邻近像素的属性值有关。两个像素点p和q之间的混合距离Dm(p,q)等于它们之间满足混合连通的通路的长度。 例3.3混合距离示例。一张3×3的二值图像模板如图38(a)所示,相似准则V={1},像素p和q的灰度值均取1。试求: 当像素s和t的灰度值如图38(b)、图38(c)、图38(d)和图38(e)取值时,像素p和q之间的混合距离Dm(p,q)。 图38混合距离示例 解: 当像素s和t的灰度值如图38(b)取值时,满足混合连通的通路如图所示,Dm(p,q)=2。 当像素s和t的灰度值如图38(c)取值时,满足混合连通的通路如图所示,Dm(p,q)=3。 当像素s和t的灰度值如图38(d)取值时,满足混合连通的通路如图所示,Dm(p,q)=3。 当像素s和t的灰度值如图38(e)取值时,满足混合连通的通路如图所示,Dm(p,q)=4。 3.3几何变换 几何变换指将图像的几何信息进行变换来获取新图像的变换方法,包括平移变换、放缩变换、旋转变换、镜像变换、剪切变换、透视变换等。二维图像和三维图像的几何变换原理基本相同,本节主要基于三维图像几何变换的原理进行介绍。 首先讨论坐标的齐次表示问题。齐次坐标是指把一个n维向量用一个n+1维向量来表示。设空间中一个点的笛卡儿坐标为(x,y,z),则其对应的齐次坐标为(Hx,Hy,Hz,H); 当H=1时,坐标(x,y,z,1)称为坐标(x,y,z)对应的规范化齐次坐标。齐次坐标的表示方式是不唯一的,引入齐次坐标的主要目的是合并矩阵运算中的乘法和加法,它提供了一种利用矩阵运算把二维、三维甚至高维空间中的一个点集从一个坐标系变换到另一个坐标系的有效方法。 几何变换的本质是矩阵运算。设空间中一个点的笛卡儿坐标为(x,y,z),基于某个条件将其变换到新的坐标(x′,y′,z′),则几何变换公式可表示为如式(317)所示的矩阵形式: P′=x′ y′ z′ 1=a11a12a13a14 a21a22a23a24 a31a32a33a34 a41a42a43a44x y z 1=AP(317) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; A称为几何变换矩阵,它唯一地确定了几何变换的结果。 图像由像素点构成,对单个点的变换可以推广到对m个点的变换,从而实现对图像的几何变换。设m个点的笛卡儿坐标如式(318)所示,该式可表示为如式(319)所示的矩阵形式,基于某个条件将其变换到如式(320)所示的新的坐标,则几何变换公式可表示为如式(321)所示的矩阵形式: P1=x1 y1 z1 1,P2=x2 y2 z2 1,…,Pm=xm ym zm 1(318) P=x1x2…xm y1y2…ym z1z2…zm 11…1(319) P′=x′1x′2…x′m y′1y′2…y′m z′1z′2…z′m 11…1(320) P′=x′1x′2…x′m y′1y′2…y′m z′1z′2…z′m 11…1=a11a12a13a14 a21a22a23a24 a31a32a33a34 a41a42a43a44x1x2…xm y1y2…ym z1z2…zm 11…1=AP(321) 式中,A称为几何变换矩阵,它唯一地确定了几何变换的结果。 3.3.1平移变换 平移变换(Translation Transformation)是一种刚体变换(RigidBody Transformation),指将图像沿某方向平移来获取新图像的变换方法。设空间中一个点的笛卡儿坐标为(x,y,z),基于平移向量(a,b,c)将其平移到新的坐标(x′,y′,z′),则平移变换公式可表示为如式(322)所示的矩阵形式: P′=x′ y′ z′ 1=100a 010b 001c 0001x y z 1=TP=x+a y+b z+c 1(322) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; a、b、c称为平移系数,它们共同构成了平移向量(a,b,c); T称为平移变换矩阵。 平移变换的效果图如图39所示。 图39平移变换 3.3.2放缩变换 放缩变换(Scale Transformation)也称为尺度变换,指将图像在某方向按比例放缩来获取新图像的变换方法。放缩变换改变了图像的尺寸,即改变了图像像素点间的距离。放缩变换一般沿坐标轴方向进行,也可分解为沿坐标轴方向进行。 设空间中一个点的笛卡儿坐标为(x,y,z),基于放缩向量(a,b,c)将其放缩到新的坐标(x′,y′,z′),则放缩变换公式可表示为如式(323)所示的矩阵形式: P′=x′ y′ z′ 1=a000 0b00 00c0 0001x y z 1=SP=ax by cz 1(323) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; a、b、c称为放缩系数,它们共同构成了放缩向量(a,b,c); 当a=b=c时,称为尺度放缩,否则,称为拉伸放缩; S称为放缩变换矩阵。 需要注意的是,当放缩系数a、b、c不为整数时,原始图像中某些像素放缩后的坐标可能不为整数,导致变换后的图像中出现“孔洞”现象,此时,需要经过取整或插值等操作来进行失真校正。未经失真校正的放缩变换效果图如图310所示,经过失真校正的放缩变换效果图如图311所示。 图310未经失真校正的放缩变换 图311经过失真校正的放缩变换 % F3_11.m I1=imread('lena.bmp'); I2=imresize(I1,0.5); I3=imresize(I1,1.5); figure,imshow(I1),xlabel('(a) 原图像'); figure,imshow(I2),xlabel('(b) 缩小后的图像'); figure,imshow(I3),xlabel('(c) 放大后的图像'); 3.3.3旋转变换 旋转变换(Rotation Transformation)是一种刚体变换(RigidBody Transformation),指将图像以某点为轴进行旋转来获取新图像的变换方法。首先考虑2D平面上一个点绕原点的旋转。设平面上一个点的笛卡儿坐标为(x0,y0),将该点 图3122D旋转变换示意图 绕原点O顺时针旋转角度a后到达新的坐标(x1,y1),则旋转变换示意图如图312所示。其中,r为该点到原点O的距离; b为旋转前r与x轴之间的夹角。 旋转前,式(324)成立: x0=rcosb y0=rsinb(324) 旋转后,式(325)成立: x1=rcos(b-a)=rcosbcosa+rsinbsina=x0cosa+y0sina y1=rsin(b-a)=rsinbcosa-rcosbsina=-x0sina+y0cosa(325) 式(325)可表示为如式(326)所示的矩阵形式: x1 y1=cosasina -sinacosax0 y0(326) 最简单的3D旋转变换是一个点绕坐标轴的旋转,包括分别绕Z、Y、X坐标轴旋转。 1. 一个点绕Z坐标轴旋转 设空间中一个点的笛卡儿坐标为(x,y,z),将该点绕Z坐标轴旋转γ角度到达新的坐标(x′,y′,z′),则旋转变换公式可表示为如式(327)所示的矩阵形式: P′=x′ y′ z′ 1=cosγsinγ00 -sinγcosγ00 0010 0001x y z 1=RγP=xcosγ+ysinγ -xsinγ+ycosγ z 1(327) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; γ为该点绕Z坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的; Rγ称为旋转变换矩阵。 2. 一个点绕Y坐标轴旋转 设空间中一个点的笛卡儿坐标为(x,y,z),将该点绕Y坐标轴旋转β角度到达新的坐标(x′,y′,z′),则旋转变换公式可表示为如式(328)所示的矩阵形式: P′=x′ y′ z′ 1=cosβ0sinβ0 0100 -sinβ0cosβ0 0001x y z 1=RβP=xcosβ+zsinβ y -xsinβ+zcosβ 1(328) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; β为该点绕Y坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的; Rβ称为旋转变换矩阵。 3. 一个点绕X坐标轴旋转 设空间中一个点的笛卡儿坐标为(x,y,z),将该点绕X坐标轴旋转α角度到达新的坐标(x′,y′,z′),则旋转变换公式可表示为如式(329)所示的矩阵形式: P′=x′ y′ z′ 1=1000 0cosαsinα0 0-sinαcosα0 0001x y z 1=RαP=x ycosα+zsinα -ysinα+zcosα 1(329) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; α为该点绕X坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的; Rα称为旋转变换矩阵。 需要注意的是,原始图像中某些像素旋转后的坐标可能不为整数,导致变换后的图像中出现“孔洞”现象,此时,需要经过取整或插值等操作来进行失真校正。未经失真校正的旋转变换效果图如图313所示,经过失真校正的旋转变换效果图如图314所示。 图313未经失真校正的旋转变换 图314经过失真校正的旋转变换 % F3_14.m I1=imread('lena.bmp'); I2=imrotate(I1,-30); I3=imrotate(I1,30); subplot(1,3,1),imshow(I1),xlabel('(a) 原始图像'); subplot(1,3,2),imshow(I2),xlabel('(b) 顺时针旋转'); subplot(1,3,3),imshow(I3),xlabel('(c) 逆时针旋转'); 一般的3D旋转变换是一个点绕任意一个中心点旋转,此时的情况相对复杂。设空间中一个点P绕一个中心点C旋转可由三个变换的级联来实现: 第一个变换将中心点C平移到坐标原点,第二个变换将点P绕坐标原点旋转,第三个变换将中心点C平移回原来的初始位置。 3.3.4镜像变换 镜像变换(Mirror Transformation)是一种刚体变换(RigidBody Transformation),包括水平镜像和垂直镜像。 1. 水平镜像 水平镜像指将图像左半部分和右半部分以图像垂直中轴线为中心进行镜像对换来获取新图像的变换方法。设空间中一个点的笛卡儿坐标为(x,y,z),将其水平镜像到新的坐标(x′,y′,z′),则水平镜像变换公式可表示为如式(330)所示的矩阵形式: P′=x′ y′ z′ 1=-100w 0100 0010 0001x y z 1=MxP=-x+w y z 1(330) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; w为图像的宽度; Mx称为水平镜像变换矩阵。 2. 垂直镜像 垂直镜像指将图像上半部分和下半部分以图像水平中轴线为中心进行镜像对换来获取新图像的变换方法。设空间中一个点的笛卡儿坐标为(x,y,z),将其垂直镜像到新的坐标(x′,y′,z′),则垂直镜像变换公式可表示为如式(331)所示的矩阵形式: P′=x′ y′ z′ 1=1000 0-10h 0010 0001x y z 1=MyP=x -y+h z 1(331) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; h为图像的高度; My称为垂直镜像变换矩阵。 镜像变换的效果图如图315所示。 图315镜像变换 % F3_15.m I = imread('lena.bmp'); I1 = flipdim(I,2); I2 = flipdim(I,1); subplot(1,3,1),imshow(I),xlabel('(a) 原始图像'); subplot(1,3,2),imshow(I1),xlabel('(b) 水平镜像图像'); subplot(1,3,3),imshow(I2),xlabel('(c) 垂直镜像图像'); 3.3.5剪切变换 剪切变换(Shear Transformation)也称为错切变换,刻画了类似四边形不稳定性的性质,包括水平剪切和垂直剪切。 1. 水平剪切 水平剪切指将图像一条水平边固定,并沿水平方向拉长图像来获取新图像的变换方法。设空间中一个点的笛卡儿坐标为(x,y,z),将其水平剪切到新的坐标(x′,y′,z′),则水平剪切变换公式可表示为如式(332)所示的矩阵形式: P′=x′ y′ z′ 1=1000 a100 0010 0001x y z 1=HxP=x ax+y z 1(332) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; a为水平剪切系数; Hx称为水平剪切变换矩阵。 2. 垂直剪切 垂直剪切指将图像一条垂直边固定,并沿垂直方向拉长图像来获取新图像的变换方法。设空间中一个点的笛卡儿坐标为(x,y,z),将其水平剪切到新的坐标(x′,y′,z′),则水平剪切变换公式可表示为如式(333)所示的矩阵形式: P′=x′ y′ z′ 1=1b00 0100 0010 0001x y z 1=HyP=x+by y z 1(333) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; b为垂直剪切系数; Hy称为垂直剪切变换矩阵。 剪切变换的效果图如图316所示。 图316剪切变换 % F3_16.m I=imread('lena.bmp'); T = [ 100 0.510 001] tform = maketform('affine',T); J1 = imtransform(I,tform); T = [10.50 010 001] tform = maketform('affine',T); J2 = imtransform(I,tform); T = [10.50 0.510 001] tform = maketform('affine',T); J3 = imtransform(I,tform); subplot(2,2,1),imshow(I),xlabel('(a) 原始图像'); subplot(2,2,2),imshow(J1),xlabel('(b) 水平剪切'); subplot(2,2,3),imshow(J2),xlabel('(c) 垂直剪切'); subplot(2,2,4),imshow(J3),xlabel('(d) 水平垂直剪切'); 3.3.6透视变换 透视变换(Perspective Transformation)指利用透视中心、像点、目标点三点共线的条件,按透视旋转定律使透视面绕透视轴旋转某一角度,破坏原有的投影光线束,仍能保持透视面上投影几何图形不变的变换方法。透视变换原理如图317所示。 图317透视变换原理 设空间中一个点的笛卡儿坐标为(x,y,z),将其透视变换到新的坐标(x′,y′,z′),则透视变换公式可表示为如式(334)所示的矩阵形式: P′=x′ y′ z′ 1=1000 0100 0010 abc1x y z 1=EP=x y z ax+by+cz+1(334) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; a、b、c为透视变换系数,透视面为ax+by+cz=0; E称为透视变换矩阵。 透视变换的效果图如图318所示。 图318透视变换 % F3_18.m I=imread('lena.bmp'); udata = [01]; vdata = [01] T = [00 10 11 01] T1 = [-42 -8-3 1-8 66] T2 = [66 1-8 -8-3 -42] tform = maketform('projective',T,T1); [B1,xdata,ydata]= imtransform(I,tform,'bicubic','udata',udata,'vdata',vdata,'size',size(I),'fill',255); tform = maketform('projective',T,T2); [B2,xdata,ydata]= imtransform(I,tform,'bicubic','udata',udata,'vdata',vdata,'size',size(I),'fill',255); subplot(1,3,1),imshow(I),axis on,xlabel('(a) 原始图像'); subplot(1,3,2),imshow(B1),axis on,xlabel('(b) 透视图像1'); subplot(1,3,3),imshow(B2),axis on,xlabel('(c) 透视图像2'); 3.3.7反变换 几何变换的反变换是指执行与对应几何变换相反操作的变换。许多几何变换都有对应的反变换,具体描述如下。 1. 平移变换的反变换 设空间中一个点的笛卡儿坐标为(x,y,z),基于平移向量(a,b,c)已被平移到新的坐标(x′,y′,z′),则平移变换的反变换公式可表示为如式(335)所示的矩阵形式: P=x y z 1=100-a 010-b 001-c 0001x′ y′ z′ 1=T-1P′=x′-a y′-b z′-c 1(335) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; -a、-b、-c称为反平移系数,它们共同构成了反平移向量(-a,-b,-c); T-1称为平移变换的反变换矩阵。 2. 放缩变换的反变换 设空间中一个点的笛卡儿坐标为(x,y,z),基于放缩向量(a,b,c)已被放缩到新的坐标(x′,y′,z′),则放缩变换的反变换公式可表示为如式(336)所示的矩阵形式: P=x y z 1=1a000 01b00 001c0 0001x′ y′ z′ 1=S-1P′=x′a y′b z′c 1(336) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; 1/a、1/b、1/c称为反放缩系数,它们共同构成了反放缩向量(1/a,1/b,1/c); S-1称为放缩变换的反变换矩阵。 3. 旋转变换的反变换 一个点绕Z坐标轴旋转的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),该点绕Z坐标轴旋转γ角度后已被旋转到新的坐标(x′,y′,z′),则旋转变换的反变换公式可表示为如式(337)所示的矩阵形式: P=x y z 1=cos(-γ)sin(-γ)00 -sin(-γ)cos(-γ)00 0010 0001x′ y′ z′ 1=R-1γP′=x′cosγ-y′sinγ x′sinγ+y′cosγ z′ 1(337) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; -γ为该点绕Z坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的; R-1γ称为旋转变换的反变换矩阵。 一个点绕Y坐标轴旋转的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),该点绕Y坐标轴旋转β角度后已被旋转到新的坐标(x′,y′,z′),则旋转变换的反变换公式可表示为如式(338)所示的矩阵形式: P=x y z 1=cos(-β)0sin(-β)0 0100 -sin(-β)0cos(-β)0 0001x′ y′ z′ 1=R-1βP′=x′cosβ-z′sinβ y′ x′sinβ+z′cosβ 1(338) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; -β为该点绕Y坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的; R-1β称为旋转变换的反变换矩阵。 一个点绕X坐标轴旋转的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),该点绕X坐标轴旋转α角度后已被旋转到新的坐标(x′,y′,z′),则旋转变换的反变换公式可表示为如式(339)所示的矩阵形式: P=x y z 1=1000 0cos(-α)sin(-α)0 0-sin(-α)cos(-α)0 0001x′ y′ z′ 1=R-1αP′=x′ y′cosα-z′sinα y′sinα+z′cosα 1(339) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; -α为该点绕X坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的; R-1α称为旋转变换的反变换矩阵。 4. 镜像变换的反变换 水平镜像的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),已被水平镜像到新的坐标(x′,y′,z′),则水平镜像变换的反变换公式可表示为如式(340)所示的矩阵形式: P=x y z 1=-100w 0100 0010 0001x′ y′ z′ 1=M-1xP′=-x′+w y′ z′ 1(340) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; w为图像的宽度; M-1x称为水平镜像变换的反变换矩阵。由此可知,水平镜像变换矩阵Mx与水平镜像变换的反变换矩阵M-1x相同。 垂直镜像的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),已被垂直镜像到新的坐标(x′,y′,z′),则垂直镜像变换的反变换公式可表示为如式(341)所示的矩阵形式: P=x y z 1=1000 0-10h 0010 0001x′ y′ z′ 1=M-1yP′=x′ -y′+h z′ 1(341) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; h为图像的高度; M-1y称为垂直镜像变换的反变换矩阵。由此可知,垂直镜像变换矩阵My与垂直镜像变换的反变换矩阵M-1y相同。 5. 剪切变换的反变换 水平剪切的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),已被水平剪切到新的坐标(x′,y′,z′),则水平剪切变换的反变换公式可表示为如式(342)所示的矩阵形式: P=x y z 1=1000 -a100 0010 0001x′ y′ z′ 1=H-1xP′=x′ -ax′+y′ z′ 1(342) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; -a为反水平剪切系数; Hx称为水平剪切变换的反变换矩阵。 垂直剪切的反变换可定义为: 设空间中一个点的笛卡儿坐标为(x,y,z),已被垂直剪切到新的坐标(x′,y′,z′),则垂直剪切变换的反变换公式可表示为如式(343)所示的矩阵形式: P=x y z 1=1-b00 0100 0010 0001x′ y′ z′ 1=H-1yP′=x′-by′ y′ z′ 1(343) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; -b为反垂直剪切系数; Hy称为垂直剪切变换的反变换矩阵。 3.3.8复合变换 复合变换是指多个几何变换复合而成的变换,其变换矩阵为各个几何变换的变换矩阵相乘,如式(344)所示。 P′=A1A2…AnP=AP(344) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; A1、A2,…,An为各个几何变换的变换矩阵; A称为复合变换矩阵。 例如,设空间中一个点的笛卡儿坐标为(x,y,z),依次进行平移、放缩、绕Z轴旋转将其变换到新的坐标(x′,y′,z′),则复合变换公式可表示为如式(345)所示的矩阵形式: P′=x′ y′ z′ 1=cosγsinγ00 -sinγcosγ00 0010 0001a2000 0b200 00c20 0001100a1 010b1 001c1 0001x y z 1 =Rγ(S(TP))=RγSTP=AP(345) 式中,P=(x,y,z,1)T为点(x,y,z)对应的规范化齐次坐标; P′=(x′,y′,z′,1)T为变换后点(x′,y′,z′)对应的规范化齐次坐标; a1、b1、c1称为平移系数,它们共同构成了平移向量(a1,b1,c1),T称为平移变换矩阵; a2、b2、c2称为放缩系数,它们共同构成了放缩向量(a2,b2,c2),S称为放缩变换矩阵; γ为该点绕Z坐标轴旋转的角度,定义为在右手坐标系下从旋转轴正向看原点是顺时针的,Rγ称为旋转变换矩阵; A称为复合变换矩阵。 例3.4平面上一个点的坐标为(1,1),基于平移向量(1,-1)依次进行平移、顺时针旋转45°和反平移变换,试求该点变换后的坐标。 解: 该变换是一个复合变换,求解过程为 P′=x′ y′ z′ 1=T-1RγTP=100-a 010-b 001-c 0001cosγsinγ00 -sinγcosγ00 0010 0001100a 010b 001c 0001x y z 1 =cosγsinγ0-a -sinγcosγ0-b 001-c 0001100a 010b 001c 0001x y z 1=cosγsinγ0acosγ+bsinγ-a -sinγcosγ0-asinγ+bcosγ-b 0010 0001x y z 1 =cos45sin450cos45-sin45-1 -sin45cos450-sin45-cos45+1 0010 00011 1 0 1=22220-1 -222201-2 0010 00011 1 0 1 =2-1 1-2 0 1 故坐标为(1,1)的点变换后的坐标为(2-1,1-2)。 3.4几何失真校正 几何失真是指图像在获取或显示过程中产生的畸变,也称几何畸变。例如,使用长焦镜头或使用变焦镜头的长焦端时容易产生枕形失真(Pincushion Distortion),使用广角镜头或使用变焦镜头的广角端时容易产生桶形失真(Barrel Distortion),使用广角镜头还容易产生透视失真(Perspective Distortion)等,如图319所示。 图319几何失真 几何失真校正是指将存在几何失真的畸变图像校正成为无几何失真的原始图像,通常以一张基准图像或一组基准点为基准去校正畸变图像。基准图像通常由没有畸变或畸变较小的摄像系统获得,畸变图像通常由畸变较大的摄像系统获得。 几何失真校正的校正方法包括直接校正法和间接校正法。校正步骤包括空间变换和灰度插值。 3.4.1直接校正法 直接校正法也称为向前映射法,它首先由空间变换函数h1(x,y)、h2(x,y)推导出反变换函数h′1(x′,y′)、h′2(x′,y′),然后依次计算出每一个畸变图像像素坐标(x′,y′)对应的校正坐标(x,y); 但(x,y)一般不为整数,无法直接将畸变坐标(x′,y′)处的灰度值赋值给校正坐标(x,y),而是将(x′,y′)处的灰度值分配给校正坐标(x,y)周围的四个像素,据此获得校正图像。直接校正法原理图如图320所示。 图320直接校正法原理图 直接校正法存在很多不足: (1) 畸变图像像素可能会被映射到校正图像之外,因此计算效率不高; (2) 校正图像像素的灰度值由畸变图像像素的贡献之和决定,导致分配时存在较多的寻址,特别是在高阶插值时; (3) 生成校正图像的像素分布不规则,容易出现像素挤压、疏密不均等现象,需要通过灰度插值方法生成规则的栅格图像。 3.4.2间接校正法 间接校正法也称为向后映射法,它假设经过校正的图像像素坐标在基准坐标系统中为等距网格的交叉点,从网格交叉点的坐标(x,y)出发计算出畸变图像上的对应坐标(x′,y′); 但(x′,y′)一般不为整数,无法直接将畸变坐标(x′,y′)处的灰度值赋值给校正坐标(x,y),而是将(x′,y′)周围像素的灰度值进行插值以得到校正坐标(x,y)处的灰度值,据此获得校正图像。间接校正法原理图如图321所示。 图321间接校正法原理图 间接校正法生成的校正图像由逐个像素通过一步插值得到,灰度插值效率较高,在实际应用中被广泛使用。间接校正法的步骤包括空间变换和灰度插值,如下式所示: x′=h1(x,y) y′=h2(x,y) f(x,y)=g(x′,y′)(346) 式中,(x,y)为无失真坐标系中像素的坐标; (x′,y′)为失真坐标系中像素的坐标; h1(x,y)、h2(x,y)为无失真坐标系向失真坐标系变换的空间变换函数; f(x,y)为无失真的原始图像; g(x′,y′)为存在失真的畸变图像。 3.4.3空间变换 空间变换是指对畸变图像像素坐标位置进行重新排列以恢复原始图像像素坐标位置的空间关系。空间变换的关键是确定空间变换函数,该函数可以通过先验知识得到,但实际中往往难以求得,通常需要采用后验校正的方法来确定。后验校正方法的基本思想是通过控制点(即基准图像上的正确像素点和畸变图像上的失真像素点)之间的对应关系来近似空间变换函数,一般采用如下式所示的多项式拟合的形式: x′=h1(x,y)=∑Ni=0∑N-ij=0aijxiyj y′=h2(x,y)=∑Ni=0∑N-ij=0bijxiyj(347) 式中,(x,y)为无失真坐标系中像素的坐标; (x′,y′)为失真坐标系中像素的坐标; h1(x,y)、h2(x,y)为无失真坐标系向失真坐标系变换的空间变换函数; aij、bij为拟合多项式的各项待定系数; N为拟合多项式的次数。 当N=1时,空间变换函数为线性拟合,可粗略地近似几何畸变,如式(348)所示; 当N=2时,空间变换函数为二次拟合,可较精确地近似几何畸变,如式(349)所示; 当N>3时,空间变换函数为高阶拟合。一般地,二次或高阶拟合被用来补偿由实际图像系统造成的空间失真。 x′=a0+a1x+a2y y′=b0+b1x+b2y(348) x′=a0+a1x+a2y+a3xy+a4x2+a5y2 y′=b0+b1x+b2y+b3xy+b4x2+b5y2(349) 为了求解方程组,可根据待定系数的个数确定基准图像和畸变图像上控制点的对数,这些控制点构成了图像中的一个多边形区域; 将图像分成一系列覆盖全图的多边形区域的集合,对每个区域都找到足够的控制点,便可计算出待拟合多项式的系数,确定空间变换函数,从而复原原始图像。失真图像和校正图像四边形区域对应点示意图如图322所示。 图322失真图像和校正图像四边形区域对应点示意图 例3.5确定线性拟合时的空间变换函数。 解: 采用线性拟合时,确定空间变换函数就是要确定a0、a1、a2、b0、b1、b2共6个多项式系数,因此需要建立6个方程联立求解。为此,找出基准图像和畸变图像上相对应的3对控制点(x1,y1)、(x2,y2)、(x3,y3)和(x′1,y′1)、(x′2,y′2)、(x′3,y′3),分别代入式(348)后得到式(350)~式(352): x′1=a0+a1x1+a2y1 y′1=b0+b1x1+b2y1(350) x′2=a0+a1x2+a2y2 y′2=b0+b1x2+b2y2(351) x′3=a0+a1x3+a2y3 y′3=b0+b1x3+b2y3(352) 写成矩阵形式后得到式(353)和式(354): x′1=a0+a1x1+a2y1 x′2=a0+a1x2+a2y2 x′3=a0+a1x3+a2y3x′1 x′2 x′3=1x1y1 1x2y2 1x3y3a0 a1 a2(353) y′1=b0+b1x1+b2y1 y′2=b0+b1x2+b2y2 y′3=b0+b1x3+b2y3y′1 y′2 y′3=1x1y1 1x2y2 1x3y3b0 b1 b2(354) 通过联立方程组求解或矩阵求逆运算,解出a0、a1、a2、b0、b1、b2共6个多项式系数,即可确定空间变换函数h1(x,y)和h2(x,y)。空间变换函数一旦确定,即可复原原始图像中此三点连线所包围的三角形区域内的各像素点; 以此类推,每三对控制点一组重复进行上述操作,即可实现全图像的空间校正。 例3.6确定二次拟合时的空间变换函数。 解: 采用二次拟合时,确定空间变换函数就是要确定a0、a1、a2、a3、a4、a5、b0、b1、b2、b3、b4、b5共12个多项式系数,因此需要建立12个方程联立求解。为此,找出原始图像和畸变图像上相对应的6对控制点(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)、(x5,y5)、(x6,y6)和(x′1,y′1)、(x′2,y′2)、(x′3,y′3)、(x′4,y′4)、(x′5,y′5)、(x′6,y′6),分别代入式(349)后得到式(355)~式(360): x′1=a0+a1x1+a2y1+a3x1y1+a4x21+a5y21 y′1=b0+b1x1+b2y1+b3x1y1+b4x21+b5y21(355) x′2=a0+a1x2+a2y2+a3x2y2+a4x22+a5y22 y′2=b0+b1x2+b2y2+b3x2y2+b4x22+b5y22(356) x′3=a0+a1x3+a2y3+a3x3y3+a4x23+a5y23 y′3=b0+b1x3+b2y3+b3x3y3+b4x23+b5y23(357) x′4=a0+a1x4+a2y4+a3x4y4+a4x24+a5y24 y′4=b0+b1x4+b2y4+b3x4y4+b4x24+b5y24(358) x′5=a0+a1x5+a2y5+a3x5y5+a4x25+a5y25 y′5=b0+b1x5+b2y5+b3x5y5+b4x25+b5y25(359) x′6=a0+a1x6+a2y6+a3x6y6+a4x26+a5y26 y′6=b0+b1x6+b2y6+b3x6y6+b4x26+b5y26(360) 合并坐标方程后得到式(361)和式(362): x′1=a0+a1x1+a2y1+a3x1y1+a4x21+a5y21 x′2=a0+a1x2+a2y2+a3x2y2+a4x22+a5y22 x′3=a0+a1x3+a2y3+a3x3y3+a4x23+a5y23 x′4=a0+a1x4+a2y4+a3x4y4+a4x24+a5y24 x′5=a0+a1x5+a2y5+a3x5y5+a4x25+a5y25 x′6=a0+a1x6+a2y6+a3x6y6+a4x26+a5y26(361) y′1=b0+b1x1+b2y1+b3x1y1+b4x21+b5y21 y′2=b0+b1x2+b2y2+b3x2y2+b4x22+b5y22 y′3=b0+b1x3+b2y3+b3x3y3+b4x23+b5y23 y′4=b0+b1x4+b2y4+b3x4y4+b4x24+b5y24 y′5=b0+b1x5+b2y5+b3x5y5+b4x25+b5y25 y′6=b0+b1x6+b2y6+b3x6y6+b4x26+b5y26(362) 写成矩阵形式后得到式(363)和式(364): x′1 x′2 x′3 x′4 x′5 x′6=1x1y1x1y1x21y21 1x2y2x2y2x22y22 1x3y3x3y3x23y23 1x4y4x4y4x24y24 1x5y5x5y5x25y25 1x6y6x6y6x26y26a0 a1 a2 a3 a4 a5(363) y′1 y′2 y′3 y′4 y′5 y′6=1x1y1x1y1x21y21 1x2y2x2y2x22y22 1x3y3x3y3x23y23 1x4y4x4y4x24y24 1x5y5x5y5x25y25 1x6y6x6y6x26y26b0 b1 b2 b3 b4 b5(364) 通过联立方程组求解或矩阵求逆运算,解出a0、a1、a2、a3、a4、a5、b0、b1、b2、b3、b4、b5共12个多项式系数,即可确定空间变换函数h1(x,y)和h2(x,y)。空间变换函数一旦确定,即可复原原始图像中此六点连线所包围的六边形区域内的各像素点; 以此类推,每6对像素点一组重复进行上述操作,即可实现全图像的空间校正。 3.4.4灰度插值 灰度插值是指对图像映射位置及其周围像素的灰度值进行插值操作以复原原始图像像素的灰度值。基于间接矫正法复原原始图像时,对原始图像中的每一个像素点(x,y),计算出畸变图像上的对应坐标(x′,y′),基于(x′,y′)处的灰度值确定原始图像中像素点(x,y)的灰度值。如果计算出的畸变图像上的对应坐标(x′,y′)为整数,则原始图像对应像素点(x,y)的灰度值与其保持一致; 如果不为整数,则需要进行灰度插值操作。常用的灰度插值方法包括最近邻插值法、双线性插值法和双三次插值法等。 1. 最近邻插值法 最近邻插值法(Nearest Neighbor Interpolation)也称为零阶插值,是指将距离映射位置最近的像素点的灰度值作为插值结果的方法,其原理图如图323所示。 图323最近邻插值法原理图 设原始图像f(x,y)中某像素点(x,y)经过变换后在畸变图像g(x′,y′)上的映射位置为(x′,y′),则最近邻插值法如下式所示: f(x,y)=g(x′,y′)=g(u′,v′)(365) 式中,像素点(u′,v′)为距离映射位置(x′,y′)最近的像素点,即u′、v′满足下式条件: x′-12<u′<x′+12 y′-12<v′<y′+12(366) 最近邻插值法简单易行,计算量小,但细微结构粗糙,容易出现块状效应。最近邻插值法的效果图如图324所示。 图324最近邻插值法效果图 % F3_24.m I=imread('lena.bmp'); subplot(2,2,1),imshow(I),xlabel('(a) 原始图像'); times=8; f=imresize(I,1/times); subplot(2,2,2),imshow(f),xlabel('(b) 失真图像'); % 方法1 g1=interp2(double(f),'nearest'); % g1=im2uint8(mat2gray(g1)); subplot(2,2,3),imshow(uint8(g1)),xlabel('(c) 最近邻插值(方法1)'); % 方法2 [m,n]=size(f); [x,y]=meshgrid(1:m,1:n); [xi,yi]=meshgrid(1:m*times,1:n*times); g2=interp2(x,y,double(f),xi/times,yi/times,'nearest'); % g2=im2uint8(mat2gray(g2)); subplot(2,2,4),imshow(uint8(g2)),xlabel('(d) 最近邻插值(方法2)'); 2. 双线性插值法 双线性插值法(Bilinear Interpolation)是指将映射位置周围4个像素点的灰度值在水平和垂直两个方向上进行插值以获取插值结果的方法,其原理图如图325所示。 图325双线性插值法原理图 设原始图像f(x,y)中某像素点(x,y)经过变换后在畸变图像g(x′,y′)上的映射位置为G(x′,y′),该映射位置周围4个像素点的坐标为A([x′],[y′])、B([x′],[y′]+1)、C([x′]+1,[y′]+1)和D([x′]+1,[y′]),符号[]表示取整,则双线性插值法计算图如图326所示。 图326双线性插值法计算图 由图可知下式成立: a=x′-[x′] b=y′-[y′](367) 因此,点E处的灰度值g(E)可表示为 g(E)=g(A)+b(g(B)-g(A))=(1-b)g(A)+bg(B)(368) 点F处的灰度值g(F)可表示为 g(F)=g(C)+b(g(D)-g(C))=(1-b)g(C)+bg(D)(369) 点G处的灰度值g(G)可表示为 g(G)=g(E)+a(g(F)-g(E))=(1-a)g(E)+ag(F) =(1-a)[(1-b)g(A)+bg(B)]+a[(1-b)g(C)+bg(D)] =(1-a)(1-b)g(A)+(1-a)bg(B)+a(1-b)g(C)+abg(D)(370) 因此,双线性插值法可表示为 f(x,y)=g(x′,y′)=g(G)(371) 当x=[x′]时,有a=0,此时,双线性插值法可表示为 f(x,y)=g(x′,y′)=g(G)=(1-b)g(A)+bg(B) =(1-y′+[y′])g([x′],[y′])+(y′-[y′])g([x′],[y′]+1)(372) 当y=[y′]时,有b=0,此时,双线性插值法可表示为 f(x,y)=g(x′,y′)=g(G)=(1-a)g(A)+ag(C) =(1-x′+[x′])g([x′],[y′])+(x′-[x′])g([x′]+1,[y′])(373) 当x=[x′],y=[y′]时,有a=0,b=0,此时,双线性插值法可表示为 f(x,y)=g(x′,y′)=g(G)=g(A)=g([x′],[y′])(374) 双线性插值法具有低通滤波器的性质,它削弱了高频信息,导致图像轮廓模糊。双线性插值法的效果图如图327所示。 图327双线性插值法效果图 % F3_27.m I=imread('lena.bmp'); subplot(2,2,1),imshow(I),xlabel('(a) 原始图像'); times=4; f=imresize(I,1/times); subplot(2,2,2),imshow(f),xlabel('(b) 失真图像'); % 方法1 g1=interp2(double(f),'linear'); % g1=im2uint8(mat2gray(g1)); subplot(2,2,3),imshow(uint8(g1)),xlabel('(c) 双线性插值(方法1)'); % 方法2 [m,n]=size(f); [x,y]=meshgrid(1:m,1:n); [xi,yi]=meshgrid(1:m*times,1:n*times); g2=interp2(x,y,double(f),xi/times,yi/times,'linear'); % g2=im2uint8(mat2gray(g2)); subplot(2,2,4),imshow(uint8(g2)),xlabel('(d) 双线性插值(方法2)'); 3. 双三次插值法 双三次插值法(Bicubic Interpolation)是指将映射位置周围16个像素点的灰度值在水平和垂直两个方向上进行插值以获取插值结果的方法,其原理图如图328所示: 图328双三次插值法原理图 设原始图像f(x,y)中某像素点(x,y)经过变换后在畸变图像g(x′,y′)上的映射位置为G(x′,y′),该映射位置周围16个像素点分别为G0,G1,…,G15,符号[]表示取整,则双三次插值法计算图如图329所示。 由图可知下式成立: a=x′-[x′] b=y′-[y′](375) 双三次插值法首先基于式(376)计算出点A0、A1、A2、A3处的灰度值,再基于式(377)计算出最终的插值结果。 g(Ai)=∑3j=0g(G4i+j)s(dy(G4i+j,G)),i=0,1,2,3(376) g(x′,y′)=∑3i=0g(Ai)s(dx(Ai,G))(377) 式中,两点P1(x1,y1)和P2(x2,y2)之间的距离函数d()定义为如式(378)所示; 理论上为了逼近最佳插值函数y=sinx/x,三次多项式函数s(x)定义如式(379)所示,其函数曲线如图330所示。 dx(P1,P2)=|x1-x2| dy(P1,P2)=|y1-y2|(378) s(x)=1-2|x|2+|x|3,|x|<1 4-8|x|+5|x|2-|x|3,1≤|x|<2 0,|x|>2(379) 图329双三次插值法计算图 图330s(x)函数曲线 % F3_30.m x=-6.3:0.01:6.3; y=sin(x)./x; subplot(1,2,1),plot(x,y);hold on; line([-8 8],[0 0]);hold on; x1=-1:0.01:1; s1=1-2.*abs(x1).*abs(x1)+abs(x1).*abs(x1).*abs(x1); subplot(1,2,1),plot(x1,s1);hold on; x2=-2:0.01:-1; s2=4-8.*abs(x2)+5.*abs(x2).*abs(x2)-abs(x2).*abs(x2).*abs(x2); subplot(1,2,1),plot(x2,s2);hold on; x3=1:0.01:2; s3=4-8.*abs(x3)+5.*abs(x3).*abs(x3)-abs(x3).*abs(x3).*abs(x3); subplot(1,2,1),plot(x3,s3),xlabel('(a) y=sin(x)/x与双三次插值函数对比');hold on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% x1=-1:0.01:1; s1=1-2.*abs(x1).*abs(x1)+abs(x1).*abs(x1).*abs(x1); subplot(1,2,2),plot(x1,s1,'k');hold on; x2=-2:0.01:-1; s2=4-8.*abs(x2)+5.*abs(x2).*abs(x2)-abs(x2).*abs(x2).*abs(x2); subplot(1,2,2),plot(x2,s2,'k');hold on; x3=1:0.01:2; s3=4-8.*abs(x3)+5.*abs(x3).*abs(x3)-abs(x3).*abs(x3).*abs(x3); subplot(1,2,2),plot(x3,s3,'k'),xlabel('(b) 双三次插值函数');hold on; plot([-2.5 2.5],[0 0],'k');plot([0 0],[-0.2 1.1],'k'); 由此可知,点A0、A1、A2、A3处的灰度值由式(380)~式(383)计算,它们可表示为如式(384)所示的矩阵形式。 g(A0)=∑3j=0g(G0+j)s(dy(G0+j,G)) =g(G0)s(dy(G0,G))+g(G1)s(dy(G1,G))+g(G2)s(dy(G2,G))+ g(G3)s(dy(G3,G)) =g(G0)s(1+b)+g(G1)s(b)+g(G2)s(1-b)+g(G3)s(2-b)(380) g(A1)=∑3j=0g(G4+j)s(dy(G4+j,G)) =g(G4)s(dy(G4,G))+g(G5)s(dy(G5,G))+g(G6)s(dy(G6,G))+ g(G7)s(dy(G7,G)) =g(G4)s(1+b)+g(G5)s(b)+g(G6)s(1-b)+g(G7)s(2-b)(381) g(A2)=∑3j=0g(G8+j)s(dy(G8+j,G)) =g(G8)s(dy(G8,G))+g(G9)s(dy(G9,G))+g(G10)s(dy(G10,G))+ g(G11)s(dy(G11,G)) =g(G8)s(1+b)+g(G9)s(b)+g(G10)s(1-b)+g(G11)s(2-b)(382) g(A3)=∑3j=0g(G12+j)s(dy(G12+j,G)) =g(G12)s(dy(G12,G))+g(G13)s(dy(G13,G))+g(G14)s(dy(G14,G))+ g(G15)s(dy(G15,G)) =g(G12)s(1+b)+g(G13)s(b)+g(G14)s(1-b)+g(G15)s(2-b)(383) g(A0) g(A1) g(A2) g(A3)=g(G0)g(G1)g(G2)g(G3) g(G4)g(G5)g(G6)g(G7) g(G8)g(G9)g(G10)g(G11) g(G12)g(G13)g(G14)g(G15)s(1+b) s(b) s(1-b) s(2-b)(384) 最终的插值结果由式(385)计算,它们可表示为如式(386)所示的矩阵形式。 g(x′,y′)=∑3i=0g(Ai)s(dx(Ai,G)) =g(A0)s(dx(A0,G))+g(A1)s(dx(A1,G))+g(A2)s(dx(A2,G))+ g(A3)s(dx(A3,G)) =g(A0)s(1+a)+g(A1)s(a)+g(A2)s(1-a)+g(A3)s(2-a)(385) g(x′,y′)=s(1+a) s(a) s(1-a) s(2-a)Tg(A0) g(A1) g(A2) g(A3)(386) 因此,双三次插值法可表示为 g(x′,y′)=s(1+a) s(a) s(1-a) s(2-a)Tg(G0)g(G1)g(G2)g(G3) g(G4)g(G5)g(G6)g(G7) g(G8)g(G9)g(G10)g(G11) g(G12)g(G13)g(G14)g(G15)s(1+b) s(b) s(1-b) s(2-b)=AGB (387) 双三次插值法计算精度高,较好地保持了图像的边缘细节,但计算量较大。双三次插值法的效果图如图331所示。 图331双三次插值法效果图 % F3_31.m I=imread('lena.bmp'); subplot(2,2,1),imshow(I),xlabel('(a) 原始图像'); times=4; f=imresize(I,1/times); subplot(2,2,2),imshow(f),xlabel('(b) 失真图像'); % 方法1 g1=interp2(double(f),'cubic'); % g1=im2uint8(mat2gray(g1)); subplot(2,2,3),imshow(uint8(g1)),xlabel('(c) 双三次插值(方法1)'); % 方法2 [m,n]=size(f); [x,y]=meshgrid(1:m,1:n); [xi,yi]=meshgrid(1:m*times,1:n*times); g2=interp2(x,y,double(f),xi/times,yi/times,'cubic'); % g2=im2uint8(mat2gray(g2)); subplot(2,2,4),imshow(uint8(g2)),xlabel('(d) 双三次插值(方法2)'); 习题 31试简述邻域的定义和类型。 32试简述邻接的定义和类型。 33试简述连接的定义和类型。 34两个图像子集P和Q如图332所示,相似准则V={1},试求: 图332图像子集P和Q (1) 子集P和子集Q是否: ①4邻接; ②8邻接。 (2) 子集P和子集Q是否: ①4连接; ②8连接; ③m连接。 (3) 子集P和子集Q是否: ①4连通; ②8连通; ③m连通。 (4) 若将子集P和子集Q以外的所有像素看成另一个子集R,试指出子集P和子集Q是否与子集R: ①4连接; ②8连接; ③m连接。 35图像子集如图333所示,试求: (1) 若相似准则V={0,1},试计算p和q之间4通路、8通路和m通路的长度; (2) 若相似准则V={1,2},试计算p和q之间4通路、8通路和m通路的长度。 36试简述欧氏距离、城区距离和棋盘距离的几何意义。 37像素之间的位置关系如图334所示,试求: 图333图像子集 图334像素位置关系 (1) 像素p和像素r之间的欧氏距离DE、城区距离D4和棋盘距离D8; (2) 像素p和像素s之间的欧氏距离DE、城区距离D4和棋盘距离D8; (3) 像素p和像素t之间的欧氏距离DE、城区距离D4和棋盘距离D8。 38像素之间的位置关系如图335所示,相似准则V ={1},试求: 图335像素位置关系 (1) 像素p和像素r之间的混合距离Dm; (2) 像素p和像素s之间的混合距离Dm; (3) 像素p和像素t之间的混合距离Dm。 39试简述几何变换的分类及含义。 310若平移变换矩阵T和尺度变换矩阵S如图336所示,试计算: 图336矩阵 (1) 对空间点P(1,2,3)先进行平移变换、再进行尺度变换所得到的结果; (2) 对空间点P(1,2,3)先进行尺度变换、再进行平移变换所得到的结果。 311若平移变换矩阵T和尺度变换矩阵S如图337所示,试计算: 图337矩阵 (1) 对空间点P(4,5,6)先进行平移变换、再进行尺度变换、最后进行反平移变换所得到的结果; (2) 对空间点P(4,5,6)先进行尺度变换、再进行平移变换、最后进行反尺度变换所得到的结果。 312试求能将图像顺时针旋转45°的旋转变换矩阵,并计算利用该矩阵对图像点P(1,0)进行旋转变换后所得到的结果。 313试简述几何失真校正的直接校正法的原理及不足。 314试简述几何失真校正的间接校正法的原理及步骤。 315试简述空间变换的概念及原理。 316试简述灰度插值的概念及基于间接校正法进行灰度插值的原理。 317试简述最近邻插值法、双线性插值法和双三次插值法的概念。