sba一个通用的稀疏光束法平差的软件包解析.docx
《sba一个通用的稀疏光束法平差的软件包解析.docx》由会员分享，可在线阅读，更多相关《sba一个通用的稀疏光束法平差的软件包解析.docx（15页珍藏版）》请在冰豆网上搜索。
sba一个通用的稀疏光束法平差的软件包解析
如果你来到这个页面来寻找一个通用的Levenberg-Marquardt算法的C/C++实现，请看levmar
引言：
本页面是关于sba，一个通用的稀疏光束法平差的C/C++软件包。
它基于GNU通用公共许可证GPL分发的。
光束法平差（BA）是作为每个基于特征的多视重建视觉算法的最后一步，用来获得最佳的三维结构和运动（如相机矩阵）参数估计。
提供初始估计，BA同时精化运动和结构参数，通过最小化观测和预测的图像点之间的投影误差。
最小化一般通过Levenberg-Marquardt（LM）算法来辅助完成。
然而，由于许多未知的因素作用于最小投影误差，一个通用的LM算法的实现（如MINPACK的lmder）当应用于BA背景下的定义的最小化问题时，会带来极高的计算代价。
幸运的是，在基本的法方程中不同的三维点和相机参数相互之间影响较小，呈现一种稀疏的块结构（如图）。
Sba利用这种稀疏的特性，使用LM算法的简化的稀疏变量来降低计算的复杂度。
Sba是通用的，因为它保证了用户对于相机和三维结构的描述参数的定义的完全控制。
因此，它事实上可以支持任何多视重建问题的显示和参数化。
比如任意投影相机，部分的或完全标定的相机，由固定的三维点进行外方位元素（即姿态）的估计，精化本征参数，等等。
用户要想在这类问题中使用sba，只需要提供合适的程序对这些问题和参数来计算估计的图像投影和他们的函数行列式（Jacobian）。
用来计算解析的函数行列式可以是手头的代码，或者使用支持符号微分的工具（如maple）生成的代码，或者通过自动微分技术获得的代码。
也可以使用近似的函数行列式，辅之以有限差分的方法。
另外，sba包含了检查用户提供的函数行列式的一致性的程序。
就我们的知识之所及，sba是第一个并且也是当前独一无二的的软件包，因为他能够不受版权限制以源代码形式放置在任何工程中。
作为sba的效率的一个指标，我们在这里说明，sba的单次测试已经涉及54台相机和5207三维点，产生了24609个图像投影。
相应的最小化问题依赖于15999个变量，sba使用非最优的BLAS在IntelP4@1.8GHzrunningLinux机器上大约7秒钟内解决。
如果没有BA的稀疏实现，那么这种规模的问题会变得非常棘手。
又一个光束法平差库，由德国斯图加特大学发布。
名字很怪，不知道全称是什么。
引言
程序DGAP实现了光束法平差的摄影测量方法，由HelmutSchmidandDuaneBrown发明。
它
基于图像和目标的几何关系的中心投影，使用最小二乘法。
特点
Camera-/self-/simultan标定，连同作者Brown,EbnerandGruen建议附加的参数。
两者可选的图像模型：
直接线性变换（DLT）和仿射变换。
不同测量位置和/或姿态数据（GPS支持的空三、直接地理参考）的集成。
精确的计算内外几何参数
测试附加参数的意义
计算协方差（新！
）和相关性
新：
对分格摄影机图像（framecameraimagery）的扩展的摄影测量模型。
新：
线扫描仪（linescanner）图像的直接地理参考。
新：
空三（aerialtriangulation，AT）的样例，GPS支持的空三并且直接地理参考。
新：
ADS-40线扫描仪图像的直接地理参考样例。
版权/许可
Copyright（C）2005DirkStallmann
Thisprogramisfreesoftware;youcanredistributeitand/ormodifyitunderthetermsofthe
GNUGeneralPublicLicenseaspublishedbytheFreeSoftwareFoundation;eitherversion2ofthe
License,or（atyouroption）anylaterversion.
Thisprogramisdistributedinthehopethatitwillbeuseful,butWITHOUTANYWARRANTY;
withouteventheimpliedwarrantyofMERCHANTABILITYorFITNESSFORAPARTICULAR
PURPOSE.SeetheGNUGeneralPublicLicenseformoredetails.
YoushouldhavereceivedacopyoftheGNUGeneralPublicLicensealongwiththisprogram;if
not,writetotheFreeSoftwareFoundation,Inc.,59TemplePlace-Suite330,Boston,MA
02111-1307,USA.
下载地址：
http:
//www.ifp.uni-stuttgart.de/publications/software/openbundle/index.en.html
依赖
库genlib2提供了各种子程序以及模板类。
可以有选择的使用LAPACK和BLAS。
LAPACK是一个Fortran77库，用来解数字线性代数中
最常见的问题。
它又依赖于基本的线性代数子程序BLAS。
BLAS也是一个Fortran77库，提供
优化的向量和矩阵操作。
两个库都是自由软件，并且拥有版权。
LAPACK和BLAS可在netlib上获得:
lib.org
为了在MSWindows2000/XP下编译DGAP，可以使用Cygnus的GNU-win32开发工具包
“Cygwin”（），版本1.3.6或者更高。
LAPACK/BLAS库在Cygwin下
依然可用。
安装
在UNIX或者Cygwin下编译
1、在同一目录下解压压缩包。
tarxzfgenlib2-$RELEASE.tar.gz
tarxzfdgap-$RELEASE.tar.gz
2、编译genlib2库：
cdgenlib2-$RELEASE
make
3、重命名这个目录为genlib2或者创建一个符号链接：
cd..
ln-sgenlib2-$RELEASEgenlib2
4、编译DGAP程序
cd..
ln-sgenlib2-$RELEASEgenlib2
为了安装DGAP只需要从主目录中复制程序dgap到/usr/binor/usr/local/bin.。
SparseBundleAdhustment1.5使用指南
Sparsebundleadjustment即稀疏集束调整，现在广泛应用于计算机视觉领域，基本成为最后优化的标准算法，就是在已经得到的初始摄像机参数和三维点数据的基础针对投影误差进行优化，得到使得均方投影误差最小意义下的Motion和Structure。
其算法的核心是利用Levenberg-Marquardt算法，由于视觉中问题的特殊性，造成矩阵稀疏，从而针对此特性进行求解。
这里介绍的是使用比较多的一个工具函数库，由希腊学者ManolisI.A.Lourakis和AntonisA.Argyros开发，网址为：
http:
//www.ics.forth.gr/lourakis/sba。
这里我主要翻译下其中接口函数的使用说明：
需要使用的结构体是structsba_crsm，定义如下：
structsba_crsm{
intnr,nc;//稀疏矩阵的行和列
intnnz;//非0元素数量
int*val;//非0元素存储空间，大小为nnz
int*colidx//非零元素的列下标，大小为nnz
int*rowptr;//val中开始一行的位置，大小为:
nr+1,rowptr[nr]=nnz
}
本结构体用于保存稀疏矩阵的元素
（根据其定义可以看出，其存储过程其实对于稀疏矩阵按行扫描，得到非零元素出现的行和列，然后用colidx存储列下标，rowptr存储数据中新的一行出现的位置，非零结果存在val中。
）
SparseBundleAdjustment通过函数sba_motstr_levmar_x（）执行，函数原型如下：
intsba_motstr_levmar_x（constintn,constintm,constintmcon,char*vmask,double*p,constintcnp,constintpnp,double*x,constintmnp,
void（*func）（double*p,structsba_crsm*idxij,int*wk1,int*wk2,double*hx,void*adata）,
void（*fjac）（double*p,structsba_crsm*idxij,int*wk1,int*wk2,double*hx,void*adata）,
void*adata,intitmax,intverbose,doubleopts[3],doubleinfo[10];
函数运行成功终止则返回达到最小需要的迭代次数，否则输入-1。
当前版本中假设观测向量的协方差矩阵为单位矩阵。
函数的核心是通过LU分解求解线性系统，使用LAPACK实现。
下面是用I,O分别表示输入和输出：
n:
3D点的数量。
（I）
m：
摄像机的数量（图像）（I）
mcon：
从第一个摄像机开始不需要修正参数的摄像机的数量。
当我们将世界坐标系与第一个摄像机的坐标对齐的时候，第一个摄像机矩阵保持为[I0]（I）
vmask：
点的可见性标志：
mask[（i-1）*m+j-1]=1如果点i在图像j中可见，否则为0.需要注意的是点和摄像机的表示从1，2开始而矩阵下标则是从0开始。
Vamsk的大小为n*m（I）；
p：
作为输入，初始参数向量P0=（a1,a2,…am,b1,b2,…bn）。
作为输出，估计的最小值，其大小为m*cnp+n*pnp（I/O）;
cnp:
定义一个摄像机的参数数量。
例如，欧氏摄像机参数化使用6个参数（3个旋转参数+3个平移参数）。
如果使用四元数表述旋转则参数的数量增加到7（4+3）。
完整的透视摄像机使用11个参数表示，如果包括全局的尺度因子则参数为12个（I）
pnp：
定义一个3D点需要的参数数量，对于欧氏几何为3，对于投影空间为4（I）
x：
观测向量X包括所有图像投影，顺序为（x11,x12,…,x1m,….xnm）.对于不可见的点，对应的xij消失。
最大为n*m*mnp.（I）
mnp：
定义一个图像点参数的数量（一般为2）（I）
func：
计算估计的参数向量的函数。
（I）
fjac：
在jac中评估
--[if!
vml]-->
--[endif]-->处的稀疏雅可比行列式
--[if!
vml]-->
--[endif]-->（I）
adata：
指向可能的额外数据的指针，传递给func，fjac。
主要是为了避免对于全局变量的直接应用（I）
itmax：
Levenberg-Marquardt迭代的最大次数。
（I）
verbose：
冗长程度。
0表示冷静的操作，大的值对应着增加冗长级别。
（I）
opts：
Levenberg-Marquardt算法中最小化参数选项，
--[if!
vml]-->
--[endif]-->，分别对应着初始衰减项的尺度因子和结束的终止门限（I）；
info：
关于最小化输出的信息，如果不需要可以设置为NULL。
（O）
info[0]：
--[if!
vml]-->
--[endif]-->初始参数估计的误差。
主要到info[0]除以所有的图像点观测的数量对应着初始均方投影误差；
info[1-4]：
（
--[if!
vml]-->
--[endif]-->，
--[if!
vml]-->
--[endif]-->，
--[if!
vml]-->
--[endif]-->，
--[if!
vml]-->
--[endif]-->），全部是在
--[if!
vml]-->
--[endif]-->计算得到的。
类似于info[0]，info[1]除以图像观测点的数量得到最终的均方投影误差；
info[5]：
总的迭代次数；
info[6]：
结束的原因：
1.过小的
--[if!
vml]-->
--[endif]-->
2过小的
--[if!
vml]-->
--[endif]-->
3迭代次数达到itmax
4增强法方程矩阵奇异，最小化过程应该从当前解重新开始采用增加的衰减项；
Info[7]:
function评价的总的次数；
Info[8]：
fjac评价的总的次数；
Info[9]:
增强法方程求解总的次数。
这个参数一般比总的迭代次数要大，由于在一次LM迭代中，可能需要多个参数进行尝试，每一个都需要对应的增强法方程的求解。
另外Sba还提供了两个函数sba_mot_levmar_x（）和sba_str_levmat_x（），分别只对相机参数和结构参数优化投影误差。
个别公式无法显示。
关于sba（sparsebundleadjustment）的30个常见问题（1-14）
sbaFAQ
Q1--什么是sba？
sba是一个C/C++软件包对广义稀疏光束平差，在GNU公共许可证下分发。
sba是通用的，提供关于定义涉及光束法平差的图像投影的参数选择和函数关系增强的灵活性。
Q2--什么是光束法平差？
假设给定一系列图像中观测到的一组对应点集相应的三维坐标的初始估计，以及关于每张图像的viewing参数的初始估计。
光束法平差（BA）是一个大的最优化的问题，包括同时精化三维结构和viewing参数（即相机姿态和可能的本征校准和径向畸变），为了获得一个在特定的假设下最优化的重建，考虑与观测的图像特征有关的噪声：
如果图像误差满足均值为零的正态分布，那么BA是最大似然法估计。
它的名字“bundles”（光束）源于每个三维特征聚焦于每个相机的光学中心，这些光学中心相对于结构和viewing参数进行最优化的调整。
sba使用Levenberg-Marquardt非线性最小二乘算法的常规实现来解决与BA相联系的稀疏的大规模的优化问题。
Q3--“稀疏”是什么意思？
由于不同的三维点和相机之间没有相互影响，BA过程中必须解决的线性系统（即法方程）包含许多零并且形成了一个稀疏的块结构特征。
这种结构可以通过避免存储和处理零元素来利用，从而获得可观的计算效益。
Q4--为什么通用的优化代码不能用于实现BA？
BA涉及大规模最小化问题的解决，一般涉及上千个变量。
优化算法通过迭代的进行函数的线性化，在当前估计的周围进行最小化并且获得线性系统（它的解确定了当前估计的一个增量）。
大多数优化代码（如，minpack,nl2sol/n2g等）假设这些线性系统是稠密的，即包含很多非零元素。
已知的几个这种系统的计算复杂度是O（n**3）。
调用BA涉及到大规模最小化问题的解（一般涉及上千个变量），所以很明显大多数通用的优化代码在应用BA时会引起运行效率降低和存储容量增大。
Q5--在哪里可以找到更多的关于光束法平差的信息？
关于BA的在基于视觉的重建中的应用的一个优秀的综述：
Triggs等人的BundleAdjustment:
AModernSynthesis。
更简短的说明，可参考bundleadjustmentarticleinWikipedia。
Q6--在哪里可以找到更多关于Levenberg-Marquardt算法的信息？
参考讲稿MethodsforNon-LinearLeastSquaresProblems,byK.Madsen,H.B.NielsenandO.Tingleff,TechnicalUniversityofDenmark,2004.
Q7--sba支持哪些类型的光束法平差？
Sba给它的用户对于描述相机和三维结构的参数定义的完全的控制。
因此，它可以支持不同的多视重建问题的实例，如投影重建，几何重建，本征相机参数精化，等等。
Sba也提供了程序处理前方交会和后方交会问题，其中相机姿态和场景结构分别保持不变。
Q8--谁在使用sba？
Sba对于计算机视觉、机器人、基于图像的图形学、摄影测量、测量、制图等领域的研究者和从业者是无价的。
它已经在全球许多实验室中使用，并且也是世界范围内以源代码方式提供的GNUGPL许可证而且授权商用的软件。
如果你想要知道更多关于sba的用途方面的信息，以下是我们列出的使用了sba的论文列表here。
Q9--使用sba时我需要什么？
为了能够使用所有的sba函数，必须安装LAPACK或者等价的库，检查lib.org/clapack的f2c'ed的免费版本。
据报告上述站点上的预编译的MSWin库损坏了，需要使用工程文件重新构建。
你也可以尝试这些预编译的MSWinLAPACK/BLAS库。
Q10--我怎么样编译sba？
首先，你必须保证在本地安装了LAPACK。
如果你必须安装它，请遵循以下安装说明包括LAPACK的分发。
第二步是编译sba本身。
tar压缩包包含Unix/Linux下使用gcc的makefile文件以及MSWin下使用VisualStudio的Makefile.vc。
请阅读这些文件中的注释以获得更多信息。
基于提供的makefile文件，它可以使用任何符合ANSI的编译器来直接编译。
Q11--为什么我得到一个未确定的外部符号dgesdd的链接错误？
可能是因为你的LAPACK是旧的版本。
Dgesdd在LAPACK3.0中引入。
如果你不想升级，你可以选择使用稍微慢一点的dgesvd。
Sba_lapack.c对于注释中的dgesvd有类似的调用。
Q12--不同的事物参数在哪里有详细的解释？
参考ICS/FORTHTR-340:
TheDesignandImplementationofaGenericSparseBundleAdjustmentSoftwarePackageBasedontheLevenberg-MarquardtAlgorithm,byM.I.A.LourakisandA.A.Argyros,2004.源代码中也包含了每个函数的每个参数的注释。
Q13--怎样把sba改编为我自己的BA变体？
如上所述，sba在选择描述相机、三维结构和图像投影的函数关系和参数时非常灵活。
因此，它可以支持许多种类的BA，包括任意投影或者仿射相机，部分地活着完全的本征标定相机，外方位元素（即姿态）估计，从特定的三维点、本征标定图像的三维重建，本征标定参数的精化，等等。
为了把sba改编为一个特定的问题，用户必须对相机和三维结构选择一个合适的参数化方法，然后提供相应的投影函数以及它的函数行列式的实现代码。
包括sba的演示程序更详细的阐述了处理几何BA问题的细节。
Q14--如何计算sba重建的初始点？
由于BA涉及一个迭代最小化问题的解，不同参数的初始估计应该提供给sba。
这种估计定义了一个初始的三维重建并且可以使用任何结构或者运动估计视觉算法来计算，如HartleyandZisserman's的书里面描述的方法。
本文全部译自http:
//www.ics.forth.gr/~lourakis/sba/faq.html
Q15--sba对outlying数据稳健吗？
短的回答：
No。
长的回答：
假设图像点特征的局部误差时零平均的正态分布，BA是最大似然法估计。
正是由于这个性质，使得BA如此强大，同时精化三维结构和viewing参数。
然而，注意，前述的假设不包括有grossoutliers（即不匹配的特征）的情况。
这些outlying特征应该在BA之前使用几何方法检测和消除。
这种方法的更多的细节，在这里。
Q16--怎样避免指定事务来计算投影函数行列式？
通过传入NULL作为函数行列式的对应参数，函数行列式在前向有限差分方法的辅助下进行计算。
然而，注意，我们不推荐这种选择：
为了获得最大的效率，我们建议提供一个函数解析的评估函数行列式。
Q17--我们如何通过投影函数和函数行列式来传入和使用我们自己的数据？
设想你想传入两个数组，一个是double型的，一个是integer型的。
一个快速并且有效的方法是使用全局变量。
避免使用全局变量的最简单的方法是声明一个结构体如下
structmydata{
doubledar[XXX];
intiar[YYY];
};
其中XXX和YYY表示合适的数组大小。
然后，定义一个结构体变量：
structmydatadata;
然后赋值：
data.dar[0]=7.0;
data.iar[0]=-17;
//etc
然后，调用合适sba程序，传入数据的地址作为参数，如：
ret=sba_motstr_levmar（...,proj,projac,（void*）&data,itmax,verbose,opts,info）;
你的proj和projac程序需要使用类型转换来解析提供的数据：
structmydata*dptr;
dptr=（structmydata*）adata;//adataispassedasvoid*
//supplieddatacannowbeaccessedasdptr->dar[0],etc
Q18--我如何检验用户提供的jacobi行列式？
文件sba_chkjac.c包括检查传入BA函数的函数行列式的正确性（即用户提供的投影函数的一致性）函数。
这些函数可以直接调用，通过将所有的sba_XXX_levmar（）和sba_XXX_levmar_x（）中的itmax的值指定为0。
这种情况下，这些函数立即返回0。
检查函数输出可能的梯度（即函数行列式的行）到stderr。
然而，我们应该注意，这些函数并不是100%安全可靠的，因为他们依赖点的评估，他们可能报告函数行列式是正确的。
解决这种问题的方法是试着检查其它点Point（s）处的函数行列式的正确性。
参考sba_chkjac.c查看更多细节。
同时，记住这些函数是为了检查解析函数行列式，而不是有限差分的数值近似。
Q19--为什么我获得一个LAPACK错误：
函数sba_Axb_Chol（）中错误的因式分解？
完整的消息是这样的：
LAPACKerror:
theleadingminoroforderXXisnotpositivedefinite,thefactorizationcouldnotbe
升级会员 