Tag Archives: 算法

演示视频之一 5路高清全景

近来发现过去发在youku上的几个demo视频误导了很多人,认为这就是我们现在的拼接技术。实际上这几个全景视频是很久前做的,跟我们目前的拼接质量和性能完全不在一个量级。但受制于手上缺乏全景设备,我们始终在demo更新上很不给力。今天我们先更新出一个,这是用一台有硬件缺陷的设备完成的,但拼接效果除掉少量瑕疵外还可以接受。

5路2048*1536半球形拼接,在12年的i5本上可以跑到100fps, 去除解码开销后cpu占用低于百分之十。实时拼接算法目前可以移植到移动平台上。

拼接效果背后是由优秀的全景相机标定算法支撑的。我们目前的量产标定算法可以在室内环境下以固定流程短时间内产生出理论最优的拼接参数,这是整个全景拼接技术的核心。

*****************************************************************

我们是一个自由团队,主要感兴趣的是计算机视觉方面新算法、技术的研发。目前全景这一块的技术已经成熟,但由于兴趣所限,我们更希望把接下来的时间用在新的挑战上,而将全景技术的产品化和推广应用托付给一个有想象力、有能力的团队,如果你/你们是这样的人,请联系我们。

QQ: 397692433

手机: 15165701250

全景摄像机量产标定算法

全景视频近期逐渐普遍化。Nokia发布了其天价的OZO全景摄像机;tencent、youku、sina这些网站也慢慢有了全景视频的身影;安防方面,巨鳄海康已经在深圳安防展展出了其研发中的全景安防摄像机,前端后端拼接全部都有,且传闻业内其他大佬也有动作。这给了我们不小的压力。

我们近来几个月闭关做全景技术的产品化,目前已经基本收官。这段时间就可以陆续看到由我们的全景算法支持的产品出现在市场上。

量产和之前作坊式的手动标定全景相机要求有着不小的差距,我们特意设计了新的标定算法来解决短时间高质量且对场地要求小的标定算法。传统的标定算法具有相当大缺陷以至于难以用到量产上,主要体现在要求使用远景,输出结果质量不稳定且难以达到理论最优结果,要求相机具有足够大的重合视野等问题,这些问题在量产标定算法中都不存在,它可以在室内场景下以固定的流程产生出理论最优的拼接效果。标定算法的量产要求应该是一些致力于全景技术团队的折戟之地,但阻碍不了我们:-P

由于精力有限,我们目前阶段只和有优秀硬件技术团队的厂商合作做全景摄像机产品。所谓优秀,是指年轻有想象力,且有让想象力落地的实力。敢而且能。当然,全景摄像机算不得新颖也算不上高门槛,在全景相机之后,我们会做一些有趣些的东西。

如果你是这样的人或团队,请联系我们。

QQ: 397692433

手机: 15165701250

Levenberg-Marquardt算法

匆匆一更。还没写完,对不起我太懒了T_T

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明,否则将追究法律责任。

Levenberg-Marquardt算法(下文简称LM算法)通常用于非线性最小二乘法的目标函数极小化。这是一个置信域方法(Trust-Region Method),为了防止步长太大而跳到非预期的局部极小值,这类算法自适应的调整步长。

假设f(x)为一个向量函数,每个分量都刻画了一个样本和模型预测之间的某种偏差,则非线性最小二乘目标函数为

(1)   \begin{eqnarray*} E&=&\transp{f}f \end{eqnarray*}

通常基于梯度的数值极小化算法,都是用一点的局域特征(如导数)建立该函数变化趋势的模型,以此来推测一个最可取的下降方向和步长。取一阶导数近似是线性模型,用切平面来近似函数在该点附近的变化趋势;取到二阶导数近似则是二阶模型,用一个抛物面来近似。当然,也可以取到更多阶,模型的几何直观也会越来越复杂,越来越精确的逼近真实的函数变化趋势,但是会付出更大的计算代价。当取到无穷阶时,模型在级数收敛区间内可以做到完全精确的刻画函数,此时这就是模型就是函数的泰勒展开了。

我们把f在任一点x附近展开f(x+\delta x)=f(x)+J(x)\delta x,这样Ex附近就是

(2)   \begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2\\ &=&|f|^2+2\transp{f}J(x)\delta x + |J(x)\delta x|^2 \end{eqnarray*}

注意Jf而非E的导数,这里我们得到了E的二阶展开,即用一个抛物面来近似目标函数E在一点附近的形态。这个二阶模型显然是一个线性最小二乘问题——这正是我们所要的,把一个非线性最小二乘问题转换为一系列线性最小二乘问题来求解。

因此一个简单直接的想法就是,每一迭代步构造这样一个线性最小二乘问题,通过极小化它下降到下一个迭代点,循环直到收敛。但这个方案显然有问题。我们的二阶展开模型只能在一个小邻域上有效,当\delta x增大时,被我们忽略掉的高阶项就会越来越明显,最后让二阶模型彻底失效。这样,直接根据二阶模型计算出来的迭代点对目标函数而言可能并非我们想要的,它的函数值可能比上一点还要高(因而目标函数不降反增);也可能尽管比上一点下降了但是却跳到了另一个我们不想要的凸区域里(因而算法会收敛到离初始猜测解更远的局部极小值点上);更严重的情况,模型甚至可能是退化的,即一种极端形式的抛物面——就像一本被半卷起的书,此时我们在弯折的方向仍能找到极小值,但在其垂直方向已经无法找到极小值(或者说处处极小)。

为了解决以上问题,我们需要给模型加入一些约束,把迭代步长约束在一个合理的范围内,保证在这个范围内模型足够有效。同时,我们也希望在二阶模型退化的时候,约束能够拯救它,保证算法总能求出一个极小值。这个约束很简单,如下

(3)   \begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2 + \lambda |C(x)\delta x|^2 \end{eqnarray*}

我们加入了一个二次约束项。单独看该项,它是一个中心在x点,开口向上的抛物面,或者说一个二次势阱。这个势阱像一个箍,把\delta x箍在x点附近的小邻域中。通过调节\lambda可以增强或减弱它对\delta x的约束:当\lambda充分大时,约束项成为主要贡献项,这时模型总有一个极小值,且步长较小;当\lambda较小时,约束项几乎可以忽略,算法以大步长快速往极小值接近。

再来看矩阵C(x). 这个矩阵通常被设置为对角的。此时,

(4)   \begin{eqnarray*} |C(x)\delta x|^2 &=& \transp{\delta x}\transp{C}C\delta x\\ &=&\sum_{i} \lambda_i \delta x_i^2 \end{eqnarray*}

\transp{C}C的本征值控制着每个方向上抛物面上升的快慢。可见,如果一个方向相应的本征值被设置的较小,那么该方向的约束项上升缓慢,约束较弱,步长跨度较大;否则反之。现实中使用的模型,目标函数对于各个参数的敏感度不同,有的参数略微调整就可以导致目标函数的强烈变化,有的则反之。C(x)这个参数矩阵使得我们可以针对各个不同的参数设置不同程度的步长约束,从而能够构造出一个稳定的算法。

(未完待补)

Notes on Conjugate Gradient Method

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明,否则将追究法律责任。

这是一个关于共轭梯度法的笔记。请大家注意的是,这是个笔记,并不是一个教程,因此少不了跳跃和欠解释的地方。对CG方法了解不多的同学请移步这里

线性方程组和极小化问题

一个关于对称矩阵A的线性方程组Ax=b等价于求解如下极小值问题:

    \[ f(x) = \frac{1}{2} \transp{x} A x - \transp{b} x \]

这很容易说明,我们微分目标函数得

(1)   \begin{eqnarray*} \td f &=& \frac{1}{2}(\td\transp{x}Ax+\transp{x}A\td x)-\transp{b} \td x\\ &=& \transp{(Ax - b)} \td x \end{eqnarray*}

所以\td f=0意味着Ax=b.

x^*为问题的解,e为偏离极小值点的位移,即

(2)   \begin{eqnarray*} Ax^*&=&b\\ e&=&\Delta x=x-x^* \end{eqnarray*}

我们定义残差或负梯度r=b-Ax,容易看出,r正比于偏离梯度零点x^*的位移e:

    \[ r=-Ae \]

最速下降法和共轭梯度法简述

最速下降法:搜索方向为本轮迭代初始点的梯度方向,搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法:搜索方向为本轮迭代初始点的梯度方向(残差方向)用A-内积做正交化,扣除掉之前所有搜索方向的分量给出,搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法的理解

首先考虑一个简单情形,即A=I时。此时f退化,且x的不同分量解耦,f的极小值问题分解为各分量上的极小值问题。我们只需要在各分量方向上找到极小值,问题就得解了。具体的说,我们从任意一个初始点出发,在方向d_0上寻找极小值,然后从这个极小值出发继续在d_1寻找d_1方向的极小,以此类推最后到达全局极小值。

A为正定对称矩阵的一般情形,其实只是上述情形中的x表述在一个新坐标系下的结果(下文详细解释)。为了求解这种情形,我们考察上述方法的求解轨迹在新坐标系里的对应。

首先,A=I时的各个正交搜索方向在新坐标系里对应着一组A-正交的方向,因此我们需要设法找出这样一组正交方向,然后各自求解这些方向上的极小值。

求解特定方向上极小值的方法不变,极小点就是这个方向上梯度正交于搜索方向时的点(原因很简单,这说明搜索方向上梯度分量为零)。

最后一点,如何找到这样一组A-正交方向?

一般方法是找一组现成的完备基,进行A-正交化。问题是计算量比较大。

记迭代到第i步时,所搜索过的这样一组A-正交的搜索方向单位矢量为\{d_0,d_1,...,d_{i-1}\},它们撑起的子空间

(3)   \begin{eqnarray*} D_i=\mathrm{span}\{d_0,d_1,...,d_{i-1}\} \end{eqnarray*}

即第i步前已经搜索过的子空间。

共轭梯度法的一个关键之处在于它可以从每次迭代的r_i中快速构造出新的搜索方向d_i,而不需要对r_i进行完整的正交化手续(即不需依次扣除r_id_0,d_1d_{i-1}方向的分量)。原因在于r_i自然的和D_{i-1}A-正交,因此构造d_i只需要从r_{i}中扣除掉d_{i-1}分量即可。

因此,共轭梯度算法有两处关键,一处为通过利用A-正交关系来把搜索问题转换到一个特殊坐标系,使得各搜索方向上的极小值问题解耦;一是利用r_iD_{i-1}是A-正交的这一规律来充分简化A-正交搜索方向的构造。

接下来我们详细论证这个算法。

第一个关键点,坐标变换和子问题解耦

(4)   \begin{eqnarray*} A&=&\transp{T}T\\ y&=&T x\\ b^{\prime}&=&T^{-T} b\\ \end{eqnarray*}

其中T^{-T}表示矩阵T求逆后做转置,很拗口,不过这个不重要。

于是我们有

(5)   \begin{eqnarray*} f(x) &=& \frac{1}{2} \transp{y} y - b^{\prime T} y\\ &=& \frac{1}{2}\sum_{n}y_n^2 - b^{\prime}_n y_n \equiv \frac{1}{2}\sum_{n} P_n(y_n) \end{eqnarray*}

可见f已经被分解为一系列独立的子问题P_n(y_n)T正是将问题解耦所用的坐标变换。

在新的\{y_n\}坐标系中,整个极小值问题可以通过在一组正交方向上依次寻找极小点来解决。但求解T需要对A做分解,并非一个计算上经济的解决方法。这一个变换的真正意义在于,我们可以通过它来找到求解轨迹在\{x_n\}坐标系中的对应。

T是线性变换,因此\{y_n\}中的一组直线仍被变换为\{x_n\}中的一组直线。\{y_n\}坐标系中,我们依次沿着一组正交方向求解极小值,因此只需要找到这组方向在\{x_n\}坐标系中的对应(一组未必正交的方向),就可以等价的直接在这组方向之上寻找极小值,而不需要经过坐标变换。

考虑这组正交方向单位矢\{y_n\}

(6)   \begin{eqnarray*} &&\transp{y_i} y_j = \delta_{ij}\\ &\Leftrightarrow& \transp{(T x_i)} T x_j=\delta_{ij}\\ &\Leftrightarrow& \transp{x_i} \transp{T}T x_j = \delta_{ij}\\ &\Leftrightarrow& \transp{x_i} A x_j = \delta_{ij} \end{eqnarray*}

可见,\{y_n\}这一组正交的方向在变换到\{x_n\}坐标系中后,虽然不是\{x_n\}中通常意义的正交矢量组,却可以在A-内积\langle \bullet,\bullet \rangle_A意义下成为一组A-正交矢量。这里的A-内积,指的是以对称矩阵A定义的内积运算

(7)   \begin{eqnarray*} \langle a,b \rangle_A \equiv \transp{a} A b \end{eqnarray*}

需要注意的是,上述结论是可逆的,如果\langle x_i,x_j \rangle_A = 0,同样也可以得出,x_i,x_j对应着\{y_n\}坐标系中两个正交的方向。因此,可以直接在\{x_n\}坐标系下,寻找一组A-正交的方向,然后在这一组方向上依次极小化目标函数,就可以找到整个问题的极小值。

第二个关键点,D_i子空间结构和正交关系

算法的更新规则为

(8)   \begin{eqnarray*} e_{i+1}&=&e_i + \alpha_i d_i\\ r_{i+1}&=&-A(e_i+\alpha_i d_i)=r_i-\alpha_i A d_i\\ d_{i}&=&r_i + \sum_{j<i} \beta_{ij} d_{j}\\ \beta_{ij}&=& - \frac{\transp{r_i}d_j}{\transp{d_j}Ad_j}\\ &=& - \frac{\langle r_i,d_j \rangle }{\langle d_j,d_j \rangle _A} \end{eqnarray*}

其中\alpha_i为第i步搜索方向上的位移长度,它的大小可以通过r_{i+1}应该和d_i正交这一要求定出

(9)   \begin{eqnarray*} &&\transp{r_{i+1}}d_i=0\\ &\Rightarrow&\transp{(e_i+\alpha_i d_i)} A d_i =0\\ &\Rightarrow&\alpha_i = - \frac{\langle d_i,e_i \rangle_A }{\langle d_i,d_i \rangle_A}= \frac{\langle d_i,r_i \rangle}{\langle d_i,d_i \rangle_A} \end{eqnarray*}

算法将参数空间划分成了一个层级结构,下面我们来看看这个结构的性质。

根据上述更新规则,可以发现D_i可以由另外几个矢量组撑起。首先,由d_i的更新规则可以发现,d_ir_id_{j<i}线性组合而成,另一方面,我们有初始条件d_0=r_0,归纳而知,d_i可以由\{r_{j\leq i}\}线性组合出来。于是我们有

(10)   \begin{eqnarray*} D_i=\mathrm{span}\{r_0,r_1,...,r_{i-1}\} \end{eqnarray*}

接着我们考虑r_i的更新规则。可以看到,D_{i+1}=\mathrm{span}\{D_{i},Ad_i\}.于是我们可以从D_1=\{d_0\}D_1=\{r_0\}构造出D_i

(11)   \begin{eqnarray*} D_i&=&\mathrm{span}\{d_0,Ad_0,A^2d_0,...,A^{i-1}d_0\} = \mathrm{span}\{d_0,A D_{i-1}\}\\ &=&\mathrm{span}\{r_0,Ar_0,A^2r_0,...,A^{i-1}r_0\} = \mathrm{span}\{r_0,A D_{i-1}\} \end{eqnarray*}

接下来我们再考虑各子空间和向量之间的正交关系。首先,e_iD_i的补空间里。这点是明显的,因为每一个优化步都会优化掉e_j在相应d_j方向上的分量,因而剩余的e_i就只位于D_i的补空间里。于是我们有

(12)   \begin{eqnarray*} e_i = \sum_{j \geq i} \delta_j d_j \end{eqnarray*}

这样ed之间的A-正交关系就很明显了。由于\langle d_i,d_j \rangle_A = \delta_{ij},而e_i只包含d_j>i分量,因而

(13)   \begin{eqnarray*} \langle e_i,d_{j<i} \rangle_A =0 \end{eqnarray*}

e_i和子空间D_i也是A-正交的。这意味着r_i正交于D_i

(14)   \begin{eqnarray*} &&\transp{e_i} A d_{j<i}=0\\ &\Rightarrow&\transp{(A e_i)} d_{j<i}=0\\ &\Rightarrow&\langle r_i,d_{j<i}\rangle =0 \end{eqnarray*}

到这里就可以解释为什么我们有\langle r_i,d_{j<i-1} \rangle_A =0了:

(15)   \begin{eqnarray*} &&\langle r_i,D_i \rangle =0\\ &\Rightarrow&\langle r_i,\mathrm{span}\{d_0,A D_{i-1}\} \rangle =0\\ &\Rightarrow&\langle r_i,A D_{i-1} \rangle =0\\ &\Rightarrow&\langle r_i,D_{i-1} \rangle_A =0 \end{eqnarray*}

至此我们完成了第二个关键点的论证,且更清楚的了解了参数空间,尤其是D_i子空间的结构。