标定程序原理注解 117.docx

标定程序原理注解 117.docx

《标定程序原理注解 117.docx》由会员分享，可在线阅读，更多相关《标定程序原理注解 117.docx（47页珍藏版）》请在冰豆网上搜索。

标定程序原理注解117

标定程序原理注解

一、标定程序的总体结构

图片的读取

标定图片的检测

标定图片相关变量的确定

世界坐标的确定

求取拐点的世界坐标

求取拐点的像素坐标

求取参数的初始值

标定参数的优化

二、图片的读取

1•标定图片的检测

（1）.检测当前文件夹下是否有能进行标定的图片存在（data_calib.m函数

1）.用matlab自带的dir函数来检验当前matlab打开的文件夹下是否有标定的图

片格式并把相应的图片数量以数组的形式存储，例如

1_bmp=dir（*bmp'、

1_ras=dir（*'ras'等;

2）.把相应图片格式的数量保存在相应的变量中，如

s_bmp=size（1_bmp,1）,

s_ras=size（1_ras,1；

tot与0的大小来判定当前

3）.把所有图片的数量保存在变量s_tot中，通过判断

文件夹下是否有图片存在：

当s_tot>0,说明有；反之，则无;

4）.如果有，则用dir函数显示当前文件夹下的所有文件名；反之，则停止进行标定。

（2）.确定标定图片的名称以及格式1）.用变量名calib_name来存储标定图片的名称，用变量format_image来存储

标定图片的格式，通过matlab自带的input函数来输入要标定的图片名称及格

式;

2）.通过format_image与已经规定的几种图片格式进行比较，如果与已规定的某

种图片格式相同，则输入图片格式有效；反之则该图片格式无效；

3）.通过调用函数lcheck_directory.m—函数来检测该文件下是否存在输入所要求

标定图片名称（简称标定名）及格式的图片；检测过程如下：

1用dir函数来检测含有标定名的文件，并把其以数组的形式存储在变量名

l中；把含有该标定名的个数存储在变量N1中，把该标定名含有的字符数存储

在变量length_name中;

如果N1V0,说明输入的标定名可能是错误的；反之，把满足上述条件的文件进行进一步检测，对变量PP取1，2，……,N1进行循环检测

calib_name进

②用函数findstr把检测到的文件名l（pp）与输入的标定名

iii中（如果标

行比较，并把标定名出现在文件名中的位置表示出来保存在变量

定名为空字符则iii=1）;同样用函数findstr把检测到的文件名与输入要求的图片格式format_image进行比较，并把format_image在文件名内中的位置表示

出来保存在变量loc_ext中；把检测到的文件名中介与calib_name和

format_image之间的数字序号存储在变量string_num中；

，则说明该文件即为要标定的文

Nima_valid中；把对应的文件

中；把对应的文件中图像格式出

③如果string_num不为空字符且iii

（1）=1

件之一；把所有满足标定的文件数量存储在变量

标记为ind_validpp存储在数组变量ind_valid现的位置存储在数组变量loc_ectension中；

注意的是上述循环如果的得到的

Nima_valid=0，则需要把图像格式变成大

写再次进行循环检测。

④图像格式的大小写都检测过后:

如果Nima_valid=0，说明输入要求的图

str2num

像格式可能是错误的；反之，把所有符合条件的文件数字序号用函数

转换后存储在数组变量indices中，数字序号用函数size求出其列数后存储在数组变量string_length中；用函数min、max分别得出indices中最小最大数,把最小的数对应的图像记为第一幅图像，用max-min+1算出标定图的总数，其值存储在变量n_ima中；通过判断string_length中最大最小的数是否相等来对

N_slots以及type_numbering进行取值，若min（string_length）=max（string_length）,则令N_slots=min（string_length），type_numbering=1；否则，令N_slots=1，type_numbering=0；

⑤数组image_numbers的第一个元素赋值为min（indices），最后一个元素赋值为max（indices），两者间元素按逐渐加一的方式依次赋值。

上面整个过程可以参照下面的图1标定流程：

开始

是

是

是

否

L标定图片格式是否一

iS定名的文件名？

甲调用size函数得l元素个数记为N1

调用checkdirectory.m

调用函数findstr检测标疋名

并把其

数存在变量stot中

否

否

调用dir函数检测当前文件夹下存在的可标定图片

（即ras、bmp、tif、pgm、ppm、jpg、jpeg格式的图片

显示当前文件下

的所有文件名

调用input函数输入要

标定图片的名称及格式

输入的图片格式与

函数寻找符合条件的图片

否存在含有

没有找到要标定的图片，

输入标定名可能是错误的

巴所有含有标定名的文件名存在数组l中，

的第k个兀素进行测试，

该元素不是要找的图片

该兀素的第一位？

函数str2num检测l元素

标定名与标定格式之间的字符

否是纯数字？

是

把所有符合条件的图片总数

保存在变量Nimavalid

图1标定流程

2•标定图片相关变量的确定

（1）.读取标定的图片数据（|ma_read_calib.m）

1）

calib_nameformat_iamge等data_calib.m，若存在就调用

.调用matlab自带函数exist来判断一些变量名如是否存在，若不存在重新进行数据输入即调用check_active_imagesm函数；

2）

第一，给含有n_ima个元素的数组

1，并且调用函数find找出

.调用check_active_image.m函数的作用：

active_images进行赋值，每一个元素的值为

active_images里的非零元素的序号存在数组

变量名进行类型拓展由int变为double型,

ind_active里；第二，对存在的一些

如center_optim、est_aIphaest_dist

等;

3）.让变量i依次取1,2,……,n_ima,进行下列循环来读取图片数据：

根据上面提到的type_numbering的取值来确定变量number_ext的取值：

当type_numbering=0时，调用函数num2str把image_number（i）转换成数字；

type_numbering=1时，通过函数sprintf把数字保存在变量number_ext中；把[calib_namenumber_ext‘fOrmat_image]构成新的图像名存在变量ima_name中，

根据图像名的format_image调用相应的函数来读入图像的数据：

当

format_image=ppm时，调用loadppm.m函数；当format_image=pgm时，调用ioadpgm.m函数；当format_image=ras时，调用readras.m函数；当format_image=bmp或tif或jpg或jpeg时，贝U用matlab自带的函数imread来读入图片的数据；把读入的图片数据存储在变量Ii中；最后把令U=li;

4）.调用函数find把images_read（images_read=active_image当images_read里某个元素对应的图片没有数据读入时，令该元素为0;反之，则不变。

）中非零元素的

序号存储在数组indread内，根据indread是否是空数组来判断是否有图像数据读入：

当ind_read为空数组时，说明没有图像数据读入，此时，令no_image_file=1；反之，说明有图像数据读入，此时，令no_image_file=0；5）.调用函数size可以得到图片I_1的行数m和列数n,并且取ny=HcaI=m和nx=Wcal=n；另外，当m=480时，令small_calib_image=O反之，令small_calib_image=1；调用函数gray获得一组图像灰度化比例数组并把其赋值给map；取active_images=images_read

（2）.通过调用函数mosaic.m|把所有标定的图片在一个界面上显示

1）.如果返回的ind_read是一个非空数组，则需要调用函数mosaic.m；反之，则

不需调用函数mosaic.m；

2）.通过调用函数mosaic.m来绘制图片显示表格，其具体过程为：

①设表格的行数为n_row，列数为n_col，调用函数sqrt算出n_ima*nx/ny的开方，并通过函数floor取小于等于且最接近它的整数，把该整数赋值给n_col;同样，调用函数ceil取大于等于且最接近n_ima/n_col的整数赋值给n_row；

2通过两次嵌套循环把所有的标定图片按比例放到上面绘制的表格内，调用函数image显示绘制的图片表格，调用函数colormap使图片转化为灰度图,

调用函数title对显示图命名，调用函数set对显示图的坐标进行属性设置。

如下图2为所有变量的关系图

图2变量关系图

三•图片的标定（Click_calib.m函数）

1•世界坐标的确定

（1）.对标定过程中一些参数的设置

1）.调用函数var2fix.m来判定dX_default、dY_default以及map等变量的值是什

么类型，如果是空字符或不明确的数值（即NaN），则清除其值；反之，就保留其现有的值；

2）.调用函数input来确定需要处理图像的数量并把其值赋给变量ima_numbers,

若ima_numbers是空字符，贝U令ima_proc=12…,n_ima；反之，贝U令

ima_proc=ima_numbers

3）.根据第一幅图像，确定标定拐点窗口的大小：

调用函数input来输入wintx的值,

若wintx的值为空字符，则令wintx=wintx_default，其中，wintx_default=max（round（nx/128）,round（ny/96））;反之，则调用函数round取最接近输入值的整数赋值给wintx;确定winty值的方法同上；由wintx和winty的值算出标定拐点窗口大小（2*wintx+1）*（2*winty+1）；

4）.调用函数input输入数值来确定计算标定格子数的方法：

若输入的值为空字符时，表示自动计算标定的格子数即令manual_squares=0反之，则表示人工计算

标定的格子数即令manual_squares=1

（2）.通过调用函数dick_ima_calib.m|进行拐点的标记

1）.调用函数image加载图片，同时调用函数colormap和函数set来把图片转换为灰度图；调用函数ginput根据光标的点击来获得相应拐点的图像坐标记为（xi,yi），其中，i为正整数；

2）.调用函数cornerfinder.m进行拐点坐标的优化，优化过程为：

①根据wintx和winty的大小构造矩阵

n12

n

设wintx=m,winty=n，贝U有

e

mask

2e

m1

m

1

e

m1

m

1

e

m1

m

e

2

e

d2

m11

m

offx

offy

以上三个矩阵都是（2m+1）*（2n+1）的矩阵;

Maxlter=10,迭代循环的初始条

②分辨率取值为resolution=0.005,迭代次数

件为：

v_extra=resolution+1,compt=0;循环迭代的条件为：

norm（v_extra）>resolution且compt

3循环迭代的每次具体运算过程为：

设刚开始确定的拐点坐标记为x1,y1，

经过优化后的对应坐标记为

X1,y1，用X1

y1来表示对应x1,y1四舍五入后的

坐标，并且构造以下矩阵:

x1

x1

x1

x1

vIx

x1

x1

x1

x1

T

x1

IIx1

y1

yJI

y1

y1

0,y1

y1

vly

y1

y1

y1

y1,y1

y1

对得到的图像区域做如下处理:

SI=conv2（conv2（SI,vlx,'same）,vly,'same）表示SI先与vix卷积再与vly卷积,

卷积后取与SI大小相同的矩阵；

SI=Sl（2:

2*wintx+4,2:

2*winty+4）;[gy,gx]=gradient（SI）gx表示SI对垂直方向

的梯度，gy表示SI对水平方向的梯度；

gx=gx（2:

2*wintx+2,2:

2*winty+2）,gy=gy（2:

2*wintx+2,2:

2*winty+2）;

经过上面对矩阵的选取后，gx，gy都是（2m+1）*（2n+1）的矩阵;

接着，进行以下矩阵运算：

px=x1+offx;

py=yi+offy;

gxx=gx.*gx.*mask;

gyy=gy*gy*mask;

gxy=gx*gy*mask;

d1=gxx.*px+gxy.*pyd2=gxy.*px+gyy.*py

由矩阵的运算可以知道：

px，py,gxx，gyy，gxy，d1，d2都是（2m+1）*（2n+1）

的矩阵；以下（i,j）表示矩阵第i行第j列的元素:

2m12n1

gxxij

i1j1

2m12n1

gxyij

i1j1

2m12n1

gyyj

i1j1

bb1

bb2

dt

cb2

y1c*

bb1

b*bb2/dt

x1a*

bb2

b*bb1/dt

2m12n1

d1ij

i1j1

2m12n1

d2ij

j1

4上式是一次迭代后得到的拐点优化坐标，上述迭代完后，令

v_extra人X]，"xi，compt=compt+1，对拐点坐标继续进行迭代优

-y1y1y1如

化，直至满足停止迭代的条件为止；停止迭代后得到拐点的最终优化坐标；

3）.调用函数plot在优化后的拐点坐标处用有红色的十字架标记，同时以该优化后的坐标为中心，画出大小为（2*wintx+1）*（2*winty+1）的红色矩形作为标记；

4）

cornefinder.M

.对四个拐角的拐点都进行上述处理后，再调用函数

对由四个拐角组成的数组坐标再次进行优化；把优化后的坐标记为：

（x1,y1）、

（x2,y2）、（x3,y3）、（x4,y4），并把（x1,y1）作为原点，第一拐点指向第四拐点的方向

为X轴方向，第一拐点指向第二拐点的方向为丫轴方向，并调用函数plot以蓝

实线依次连接四个拐角点，注明原点0,X轴以及丫轴；

5）.调用函数input输入沿X轴标定格子的数目并把该数值赋值给n_sq_x:

若

n_sq_x为空字符，则说明n_sq_x=n_sq_x_default反之，贝Un_sq_x为输入的值；n_sq_y的确定方法同n_sq_x；

6）.调用函数input输入沿X轴方向每个格子的大小（单位mm）并把该数值赋值给dX:

若dX为空字符，则说明dX=dX_default;反之，则dX为输入值大小；dY的确定方法与dX相同；

2.求取拐点的像素坐标

由摄像机标定的原理可知，标定点的像素坐标与世界坐标存在如下的代数关系:

xPXw

（1）

其中，x表示标定点的像素坐标，Xw表示标定点的世界坐标，P为一个3*4的

Pl,P2,P3,P4，

矩阵；若令x（X1，y1,1）T，XwXw1,Yw1,Zw1,1T，P

其中pi为3*1的向量，i1,2,3,4；则

（1）式可以写成:

XP1*Xw1P2*Yw1P3*Zw1P4

（2）

当Zw1始终为零时，

（2）式可以写为:

xP1*Xw1p2*Yw1P4⑶

令Hp1,P2,P4,则

（1）式可写为:

xHXw

（4）

又因为H表示的是齐次坐标之间的关系，所以H的最后一个元素为1.

h11h12h13

h21h22h23，其中h331；,则⑷式可以写为:

h31h32h33

X1

h11Xw1

h12Yw1

h13

h21Xw1

h22Yw1

h23（5）

1

h31Xw1

h32Yw1

h33

若已知n个点，则可写成的矩阵形式:

Lh0,其中,

h1

h2

Xw1

Yw1

1

0

0

0

X1Xw1

X1Yw1

X1

h3

0

0

0

Xw1

Yw1

1

y1Xw1

讨仏1

y1

h4

L

h

h5

X八wn

Ywn

1

0

0

0

xX入n八wn

xY

八n■wn

Xn

h6

0

0

0

X八wn

Ywn

1

ynXwn

ynYwn

yn

h7

h8

h9

如果测量像素坐标无误差时，则Lh0，因此线性算法是求方程组Lh0的非

零解。

当测量数据存在误差时，Lh0，即方程组Lh0不可能有非零解，通

常用最小二乘方法估计模型的解，即最小二乘解（对L作奇异值分解LUDVt,

则V的最后一个列向量就是方程Lh0的最小二乘解）。

这就是直接线性方法;有时也称为最小二乘估计。

直接线性方法的主要不足之处是测量数据的较小误差可能会导致较坏的估计结果。

发生这种情况的主要原因有：

（1）测量矩阵的元素是测量数据的非线性

函数，从而放大了测量误差，使得测量矩阵严重地偏离真实矩阵，因此导致直接线性方法的估计结果与真实值有很大偏差。

（2）测量矩阵可能有较大的条件数，—

导致估计结果的不稳定性。

对上述第一个问题解决的方法是采用因子化线性方法即是先将直接线性方

法的测量矩阵分解为若干个因子矩阵的乘积（如取H（Hnorm）1*Hrem），使得

因子矩阵的元素是某一幅图像的测量数据值或者是常量在，旨在降低测量矩阵元素关于测量数据的非线性，以抑制测量误差对估计结果的影响;

对上述第二个问题解决的方法是采用归一化线性方法即对测量数据作适当的归一化变换，使得变换后的数据所对应的测量矩阵有好的条件数，以提高数值

计算的稳定性；具体地说，归一化分为以下两步：

第一步位移变换即对测量数据

作坐标位移变换，使得原始数据的重心是新坐标系的原点；第二步伸缩变换即再对数据作坐标缩放变换，使得位移变换后的数据点分布在以重心为圆心半径为

jn的球内，n是测量数据的维数。

如对xHXw变换为HnormXHnormHXw.

（1）.通过函数|click_ima_calib.m获得拐点的图像像素坐标

1）.根据输入的X轴和丫轴的格子数，构造下面矩阵（设X轴格子数n_sq_x=m1,丫

轴格子数n_sq_y=n1）:

pts0

1

m1

2

m1

3

m1

4

m1

5

m1

丄

n1

1

1

m1

丄

n1

1

2

m1

丄

n1

1

把四个顶点的坐标进行以下标记:

a00x1y1

1T，a10

1T

a11x3y3

1T，a01

y2

1T

2）.调用函数[Homo,Hnorm,inv_Hnorm]=compute_homography（[aOOa10a11

a01],[0110;0011;1111]）

就可以求得矩阵Homo求解其具体过程为:

①根据输入矩阵得出下列表达式:

mxx

i，myy

14.

4i1yi

Scxx1

4j1

mxx，Scyy

14

4j1

yjmyy

②构造出矩阵Hnorn和inv」

Hnorm如下:

1

Scxx

Hnorm0

1

Scyy

0

mxx

Scxx

myy

Scyy

1

Scxx

mxx

invHnorm

Scyy

0

myy

1

③根据以上矩阵构造矩阵

L:

x4mxx

Scxx

y4myy

Scyy

x3mxx

Scxx

y3myy

Scyy

3mxx

x

Scxx

y3myy

Scyy

x2mxx

Scxx

y2myy

Scyy

x1mxx

Scxx

y1myy

Scyy

x4mxx

Scxx

y4myy

Scyy

x3mxx

Scxx

y3myy

Scyy

x2mxx

Scxx

y2myy

Scyy

对矩阵L进行奇异值分解:

[U,S,V]

=svd（L）

其中V为9*9矩阵,

把V中第九

列的元素取出来并且对第九个元素进行归一化处理,

把第九列重新塑造成3*3矩

阵记为Hrem贝U有Homo=inv_Hnorm*Hrem3）.获得拐点的图像像素坐标，根据计算公式:

XX=Homo*pts

就可以得到拐点的初始坐标，再经过前面介绍过拐点坐标优化函数

Cornerfinder.m

就可以得到最终的像数坐标即

x=cornerfinder（XX,l,winty,wintx）-1

其中，减去1是为了把原点带回到（0,0）,因为matlab的像素原点坐标为（1,1）;

2.求取拐点的像素坐标

（1）.通过函数|click_ima_calib.m|获得拐点的世界坐标

1）.获得拐点的世界坐标，根据建立的世界坐标和输入的X轴和丫轴的格子数以

及每个格子的单位长度就可以算出对应标定点的世界坐标:

0.10.2

0.1F

0.2

0.1*mrt

X0.1*n1

0.1*n10.1*n1

0.1*n10.1*n11

0.1*n11

plot在所有标定的拐点处用红色十字架标记，并以该拐点为中心画

2）.调用函数

出大小为（2*wintx+1）*（2*winty+1）的紫色矩形作为标记;

3）.把每一幅标定图的所有相关参数X、x、n_sq_x、n_sq_ywintx、winty、dX、dY都以数组的形式保存在数组变量string.save里;

四.标定参数的获取

1.求取参数的初始值

（1）.调用函数is3D.m来判定标定的图片是否是3D的

调用函数size算出图像的的拐点数记为Np,用函数mean求出所有拐点坐标X沿X轴、丫轴、Z轴各自的均值，并把该点记为X_mean,接着令

Y=X-X_mean*on