标定程序原理注解 117.docx

1、标定程序原理注解 117标定程序原理注解一、标定程序的总体结构图片的读取标定图片的检测标定图片相关变量的确定世界坐标的确定求取拐点的世界坐标求取拐点的像素坐标求取参数的初始值标定参数的优化二、图片的读取1标定图片的检测 (1).检测当前文件夹下是否有能进行标定的图片存在(data_calib.m函数1).用matlab自带的dir函数来检验当前matlab打开的文件夹下是否有标定的图片格式并把相应的图片数量以数组的形式存储，例如1_bmp=dir( *bmp 、1_ras=dir( *ras等;2).把相应图片格式的数量保存在相应的变量中，如s_bmp=size(1_b mp ,1),s_ra

2、s=size(1_ras,1；tot与0的大小来判定当前3).把所有图片的数量保存在变量 s_tot中，通过判断文件夹下是否有图片存在：当s_tot0,说明有；反之，则无;4).如果有，则用dir函数显示当前文件夹下的所有文件名；反之，则停止进行标 定。(2).确定标定图片的名称以及格式 1).用变量名calib_name来存储标定图片的名称，用变量 format_image来存储标定图片的格式，通过 matlab自带的input函数来输入要标定的图片名称及格式;2).通过format_image与已经规定的几种图片格式进行比较，如果与已规定的某种图片格式相同，则输入图片格式有效；反之则该图片

3、格式无效；3).通过调用函数lcheck_directory.m 函数来检测该文件下是否存在输入所要求标定图片名称(简称标定名)及格式的图片；检测过程如下：1用dir函数来检测含有标定名的文件，并把其以数组的形式存储在变量名l中；把含有该标定名的个数存储在变量 N1中，把该标定名含有的字符数存储在变量length_name 中;如果N1V0,说明输入的标定名可能是错误的；反之，把满足上述条件的文件 进行进一步检测，对变量PP取1，2，,N1进行循环检测calib_name 进用函数findstr把检测到的文件名l(pp)与输入的标定名iii中(如果标行比较，并把标定名出现在文件名中的位置表示出

4、来保存在变量定名为空字符则iii=1);同样用函数findstr 把检测到的文件名与输入要求的 图片格式format_image进行比较，并把format_image在文件名内中的位置表示出来保存在变量loc_ext中；把检测到的文件名中介与calib_name和format_image之间的数字序号存储在变量 string_num中；，则说明该文件即为要标定的文Nima_valid中；把对应的文件中；把对应的文件中图像格式出如果string_num 不为空字符且iii(1)=1件之一；把所有满足标定的文件数量存储在变量标记为ind_valid pp存储在数组变量ind_valid 现的位置存

5、储在数组变量loc_ecte nsion 中；注意的是上述循环如果的得到的Nima_valid=0，则需要把图像格式变成大写再次进行循环检测。图像格式的大小写都检测过后:如果 Nima_valid=0，说明输入要求的图str2 num像格式可能是错误的；反之，把所有符合条件的文件数字序号用函数转换后存储在数组变量 indices 中，数字序号用函数 size 求出其列数后存储在 数组变量stri ng_le ngth 中；用函数min、max分别得出in dices 中最小最大数, 把最小的数对应的图像记为第一幅图像，用max-min+1算出标定图的总数，其值 存储在变量 n_ima 中；通过

6、判断 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 标定流

7、程：开始是是是否L标定图片格式是否一iS定名的文件名？甲调用size函数得l元素个数记为N1调 用 check directory.m调用函数findstr检测标疋名并把其数存在变量s tot中否否调用dir函数检测当前文件夹下存在的可标定图片（即 ras、 bmp、tif、 pgm、 ppm、jpg、jpeg 格 式的图 片显示当前文件下的所有文件名调用input函数输入要标定图片的名称及格式输入的图片格式与函数寻找符合条件的图片否存在含有没有找到要标定的图片，输入 标定名可能是错误的巴所有含有标定名的文件名存在数组l中，的第k个兀素进行测试，该元素不是要找的图片该兀素的第一位？函数str2

8、num检测l元素标定名与标定格式之间的字符否是纯数字？是把所有符合条件的图片总数保存在变量Nima valid图1标定流程2标定图片相关变量的确定(1).读取标定的图片数据(|ma_read_calib.m)1)calib_name format_iamge 等 data_calib.m，若存在就调用.调用matlab自带函数exist来判断一些变量名如 是否存在，若不存在重新进行数据输入即调用 check_active_imagesm 函数；2)第一，给含有 n_ima个元素的数组1，并且调用函数 find找出.调用 check_active_image.m 函数的作用： active_im

9、ages进行赋值，每一个元素的值为active_images里的非零元素的序号存在数组变量名进行类型拓展由int变为double型,in d_active里；第二，对存在的一些如 center_optim、 est_aIpha est_dist等;3).让变量i依次取1,2,n_ima,进行下列循环来读取图片数据：根据上面提 到的type_numbering 的取值来确定变量 number_ext的取值：当 type_numbering=0 时，调用函数 num2str 把 image_number(i)转换成数字；type_numbering=1时，通过函数 sprintf把数字保存在变量

10、number_ext中；把 calib_name number_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中；最

11、后把令U=li;4).调用函数 find 把 images_read(images_read=active_image当 images_read里某个 元素对应的图片没有数据读入时， 令该元素为0;反之，则不变。)中非零元素的序号存储在数组ind read内，根据ind read是否是空数组来判断是否有图像数据 读入：当ind_read为空数组时，说明没有图像数据读入，此时，令no_image_file=1 ； 反之，说明有图像数据读入，此时，令 n o_image_file=0； 5).调用函数size可以得到图片I_1的行数m和列数n,并且取ny=HcaI=m和 nx=Wcal=n ；另外

12、，当 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 的开方，并通过函数

13、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

14、等变量的值是什么类型，如果是空字符或不明确的数值(即NaN)，则清除其值；反之，就保留其 现有的值；2).调用函数input来确定需要处理图像的数量并把其值赋给变量 ima_numbers,若 ima_numbers 是空 字符，贝U令 ima_proc=12 ,n_ima ； 反之，贝U令ima_p roc=ima_ nu mbers3).根据第一幅图像，确定标定拐点窗口的大小：调用函数input来输入wintx的值,若wintx的值为空字符，则令 wintx=wintx_default，其中， wintx_default=max(round(nx/128),round(ny/96);反之，

15、则调用函数 round取最接 近输入值的整数赋值给wintx ;确定winty值的方法同上；由wintx和winty的值算 出标定拐点窗口大小(2*wi ntx+1)*(2*wi nty+1)；4).调用函数input输入数值来确定计算标定格子数的方法：若输入的值为空字符 时，表示自动计算标定的格子数即令 manual_squares=0反之，则表示人工计算标定的格子数即令 manual_squares=1(2).通过调用函数dick_ima_calib.m |进行拐点的标记1).调用函数image加载图片，同时调用函数colormap和函数set来把图片转换为 灰度图；调用函数ginput根

16、据光标的点击来获得相应拐点的图像坐标记为 (xi,yi)， 其中，i为正整数；2).调用函数cornerfinder.m进行拐点坐标的优化，优化过程为： 根据wintx和winty的大小构造矩阵n 1 2n设 wintx=m , winty=n，贝U有emask2 em 1m1em 1m1em 1me2ed 2m 1 1moffxoffy以上三个矩阵都是(2m+1)*(2n+1)的矩阵;Maxlter=10,迭代循环的初始条分辨率取值为resolution=0.005,迭代次数件为：v_extra=resolution+1,compt=0;循环迭代的条件为：norm(v_extra)resol

17、ution 且 comp tMaxIter;3循环迭代的每次具体运算过程为：设刚开始确定的拐点坐标记为 x1,y1，经过优化后的对应坐标记为X1 ,y1，用 X1,y1来表示对应x1,y1四舍五入后的坐标，并且构造以下矩阵:x1x1x1x1vIxx1x1x1x1T,x1II x1y1yJIy1y10,y1y1vlyy1y1y1y1 ,y1y1对得到的图像区域做如下处理:SI = conv2(conv2(SI,vlx,same),vly,same)表示 SI先与vix卷积再与 vly卷积,卷积后取与SI大小相同的矩阵；SI = Sl(2:2*wintx+4,2:2*winty+4) ; gy,g

18、x = 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 .* py d2= gxy .* px + g

19、yy .* py由矩阵的运算可以知道：px，py, gxx，gyy，gxy，d1，d2都是(2m+1)*(2n+1)的矩阵；以下(i,j)表示矩阵第i行第j列的元素:2m 12n 1gxx iji 1 j 12m 12n 1gxy iji 1 j 12m 12n 1gyy ji 1 j 1bb1bb2dtc b2y1 c*bb1b*bb2/dtx1 a*bb2b* bb1 /dt2m 12 n 1d1 iji 1 j 12m 12 n 1d2 ijj 14上式是一次迭代后得到的拐点优化坐标，上述迭代完后，令v_extra 人 X ， xi，compt=compt+1，对拐点坐标继续进行迭代优-

20、 y1 y1 y1 如化，直至满足停止迭代的条件为止；停止迭代后得到拐点的最终优化坐标；3).调用函数plot在优化后的拐点坐标处用有红色的十字架标记，同时以该优化 后的坐标为中心，画出大小为(2*wintx+1)*(2*winty+1)的红色矩形作为标记；4)cornefin der.M.对四个拐角的拐点都进行上述处理后，再调用函数对由四个拐角组成的数组坐标再次进行优化； 把优化后的坐标记为：(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)，并把(x1,y1)作为原点，第一拐点指向第四拐点的方向为X轴方向，第一拐点指向第二拐点的方向为 丫轴方向，并调用函数plot以蓝实线依次连

21、接四个拐角点，注明原点 0, X轴以及丫轴；5).调用函数input输入沿X轴标定格子的数目并把该数值赋值给 n_sq_x:若n_sq_x为空字符，则说明n_sq_x=n_sq_x_default反之，贝U n_sq_x为输入的值； n_sq_y的确定方法同n_sq_x；6).调用函数input输入沿X轴方向每个格子的大小(单位mm)并把该数值赋值给 dX:若dX为空字符，则说明dX=dX_default;反之，则dX为输入值大小；dY的 确定方法与dX相同；2.求取拐点的像素坐标由摄像机标定的原理可知，标定点的像素坐标与世界坐标存在如下的代数关系:x PXw (1)其中，x表示标定点的像素坐

22、标，Xw表示标定点的世界坐标，P为一个3*4的Pl, P2,P3, P4，矩阵；若令 x (X1，y1,1)T，Xw Xw1,Yw1,Zw1,1T，P其中pi为3*1的向量，i 1,2,3,4 ；则(1)式可以写成:X P1 * Xw1 P2 * Yw1 P3 * Zw1 P4(2)当Zw1始终为零时，(2)式可以写为:x P1*Xw1 p 2 *Yw1 P4 令H p1,P2,P4 ,则(1)式可写为:x HXw(4)又因为H表示的是齐次坐标之间的关系，所以 H的最后一个元素为1.h11 h12 h13h21 h22 h23 ，其中h33 1 ；,则式可以写为:h31 h32 h33X1h1

23、1Xw1h12Yw1h13h21Xw1h22Yw1h23 (5)1h31Xw1h32Yw1h33若已知n个点，则可写成的矩阵形式:Lh 0,其中,h1h2Xw1Yw11000X1Xw1X1Yw1X1h3000Xw1Yw11y1Xw1讨仏1y1h4L,hh5X 八wnYwn1000x X 入n八wnx Y八 n wnXnh6000X 八wnYwn1ynXwnynYwnynh7h8h9如果测量像素坐标无误差时，则 Lh 0，因此线性算法是求方程组Lh 0的非零解。当测量数据存在误差时，Lh 0，即方程组Lh 0不可能有非零解，通常用最小二乘方法估计模型的解，即最小二乘解 （对L作奇异值分解L UD

24、Vt,则V的最后一个列向量就是方程 Lh 0的最小二乘解）。这就是直接线性方法; 有时也称为最小二乘估计。直接线性方法的主要不足之处是测量数据的较小误差可能会导致较坏的估 计结果。发生这种情况的主要原因有：（1）测量矩阵的元素是测量数据的非线性函数，从而放大了测量误差，使得测量矩阵严重地偏离真实矩阵， 因此导致直接 线性方法的估计结果与真实值有很大偏差。（2）测量矩阵可能有较大的条件数，导致估计结果的不稳定性。对上述第一个问题解决的方法是采用因子化线性方法即是先将直接线性方法的测量矩阵分解为若干个因子矩阵的乘积 （如取H （Hnorm） 1 * Hrem），使得因子矩阵的元素是某一幅图像的测量

25、数据值或者是常量在， 旨在降低测量矩阵元 素关于测量数据的非线性，以抑制测量误差对估计结果的影响;对上述第二个问题解决的方法是采用归一化线性方法即对测量数据作适当 的归一化变换，使得变换后的数据所对应的测量矩阵有好的条件数 ，以提高数值计算的稳定性；具体地说，归一化分为以下两步：第一步位移变换即对测量数据作坐标位移变换，使得原始数据的重心是新坐标系的原点； 第二步伸缩变换即再 对数据作坐标缩放变换，使得位移变换后的数据点分布在以重心为圆心半径为jn的球内，n是测量数据的维数。如对x HXw变换为HnormX HnormHXw.(1).通过函数|click_ima_calib.m获得拐点的图像像

26、素坐标1).根据输入的X轴和丫轴的格子数，构造下面矩阵(设X轴格子数n_sq_x=m1, 丫轴格子数n_sq_y=n1):pts 01m12m13m14m15m1丄n111m1丄n112m1丄n11把四个顶点的坐标进行以下标记:a00 x 1 y 11 T，a101Ta11 x 3 y 31 T，a01y21T2).调用函数Homo,Hnorm,inv_Hnorm = compute_homography(aOO a10 a11a01,0 1 1 0;0 0 1 1;1 1 1 1)就可以求得矩阵Homo求解其具体过程为:根据输入矩阵得出下列表达式:mxxi，myy1 4 .4i1yiScxx

27、 14 j 1mxx ，Scyy1 44j1y j myy构造出矩阵Hnorn和invHnorm如下:1ScxxHnorm 01Scyy0mxxScxxmyyScyy1Scxxmxxinv HnormScyy0myy1根据以上矩阵构造矩阵L:x 4 mxxScxxy 4 myyScyyx 3 mxxScxxy 3 myyScyy3 mxxxScxxy 3 myyScyyx 2 mxxScxxy 2 myyScyyx 1 mxxScxxy 1 myyScyyx 4 mxxScxxy 4 myyScyyx 3 mxxScxxy 3 myyScyyx 2 mxxScxxy 2 myyScyy对矩阵L

28、进行奇异值分解:U,S,V=svd(L),其中V为9*9矩阵,把V中第九列的元素取出来并且对第九个元素进行归一化处理,把第九列重新塑造成 3*3矩阵记为 Hrem 贝U有 Homo=inv_Hnorm*Hrem 3).获得拐点的图像像素坐标，根据计算公式:XX = Homo* pts就可以得到拐点的初始坐标，再经过前面介绍过拐点坐标优化函数Cornerfi nder.m就可以得到最终的像数坐标即x= cornerfi nder(XX,l,wi nty,wi ntx)-1其中，减去1是为了把原点带回到(0,0),因为matlab的像素原点坐标为(1,1);2.求取拐点的像素坐标(1).通过函数|

29、click_ima_calib.m |获得拐点的世界坐标1).获得拐点的世界坐标，根据建立的世界坐标和输入的 X轴和丫轴的格子数以及每个格子的单位长度就可以算出对应标定点的世界坐标:0.1 0.20.1F0.20.1*mrtX 0.1*n10.1*n1 0.1*n10.1*n1 0.1* n1 10.1* n1 1plot在所有标定的拐点处用红色十字架标记，并以该拐点为中心画2).调用函数出大小为(2*wintx+1)*(2*winty+1)的紫色矩形作为标记;3).把每一幅标定图的所有相关参数 X、x、n_sq_x、n_sq_y wintx、winty、dX、 dY都以数组的形式保存在数组变量 string.save里;四.标定参数的获取1.求取参数的初始值(1).调用函数is3D.m来判定标定的图片是否是3D的调用函数size算出图像的的拐点数记为 Np,用函数mean求出所有拐点坐 标X沿X轴、丫轴、Z轴各自的均值，并把该点记为 X_mean,接着令Y=X-X_mean*on