特色文章

开卖了

Shenzhen Civic Center

Shenzhen Civic Center

我们和客户合作研发的安防类多目全景网络摄像机今年上半年已经完成产品化,开卖了。

目前的产品包括4目180度,5目半球,8目半球,以及8目车载全景等,目前最高分辨率为2400W像素,30fps,很快将会支持到4000W像素。

产品支持多种全景类型(沉浸式半球/各种投影类型的条带展开),浏览时多码流实时切换,对PC端硬件要求低,最低支持到OpenGL 2.1.

欲知更多细节,欢迎来电。

1500w像素全景监控演示视频

全景视频技术安防行业应用案例

5 6 13 14

全景视频技术的产品化之路

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

甚嚣尘上的VR炒作终于在今年平静了,这大概意味着VR技术开始进入技术成熟度曲线的第三个时期:行业的公众关注进入低谷,人们开始冷静客观评估技术的适用范围和潜力,并逐步发现有效的经营模式。

技术成熟度曲线

VR时代的到来是不可避免的,或者说它已经到来,只是还没有推到大众面前。另外,真正具有想象力和冲击力的新技术乃是紧随其后的AR,这一点可能并不像公众预期的那样。这个时代需要由一系列扎实漂亮的产品撑起(不是概念,不是Demo,不是DIY,是产品),我们这次来谈谈全景摄像机的产品化之路上有哪些曲折和挑战。当然,全景摄像机本身并非仅限于VR应用,我们也要包括安防应用。

安防监控领域

泛泛来说有两种全景视频实时拼接方案,即前端(机内)拼接后端(PC/手机)拼接。在安防领域也是如此。前端拼接直接由全景摄像机输出拼接完成的全景帧,具有很好的兼容性,可以直接像一台普通IPC一样接入旧有系统;而后端拼接是将全景摄像机看做独立的多路IPC,同时接入监控PC服务器,由PC完成实时拼接和监看。后端拼接的优势在于可以完成极高分辨率(目前我们的后端方案全景监控分辨率最高已经有9600万像素)的全景监控,但兼容性不好,需要将全景拼接SDK嵌入平台软件,不能做到“即插即用”。

从实现上来说,大概有如下几种:FPGA/DSP/CUDA/OpenGL/CPU. 前两种用于前端拼接,FPGA的开发和维护都有较高代价,CUDA和OpenGL方案具有最高的处理能力,CPU方案除非无法选择否则是应该排除的。在前端拼接方案里,还要考虑编码问题,全景帧动辄数千万的分辨率编码并不是一个简单问题。这里我们主要谈我们自己比较熟悉的CUDA/OpenGL方案。

安防监控领域对于全景摄像机有一些特殊需求。对于后端拼接全景,其拼接参数应当保存在设备之中,由设备传给平台软件完成实时监看的初始化流程,而平台软件上则对实时拼接的效率,全景模块与其他设备如球机的互动都有颇多要求,我们简单罗列如下。

  1. 拼接参数应该是一个很小(几k到几十k)的文件,方便写入设备及在网络上传输;
  2. 灵活的裁剪/融合算法。安防全景细分需求繁多,催生大量不同类型的设备,不同目数,不同镜头,不同摆位都会导致非常不一样的应用场景,算法需要在任何场景下都能够完美融合,不产生图像瑕疵。
  3. 全景拼接模块应当支持实时的子码流/主码流切换。平台软件实时监看几十上百路网络摄像机,区分子码流/主码流非常重要,这样在小视图模式下采用子码流,而大视图下自动切换到主码流,既保证了性能又保证了操作体验。
  4. 全景拼接模块在子码流输入下能够同时完成多达几十路的实时拼接播放。
  5. 应该有多种全景投影模式。除了常见的球面/柱面展开,碗形交互式展开在安防监控领域颇受欢迎,如图所示。各种投影模式之间应该能够实时切换。
  6. ROI局域放大。
  7. 与球机联动,全景纵览全局,球机实现局域放大。
  8. 全景像素坐标到输入图像像素坐标的正反向投影。

碗形交互式全景

解决以上需求的任务并非平凡。简单一例,子/主码流实时切换中,除非子码流与主码流具有同样的视野,否则无法在不重新初始化算法的前提下完成切换,这要求算法具备瞬间初始化完毕的性能。同样的性能要求也出现在实时全景投影模式切换中。

从一个算法到一个成熟产品的道路是长远的。行业里很多学术型团队最终败在懂算法不懂软件工程,无法将一个Demo级的算法提升成一款结构良好,功能灵活,充分解决行业内需求的算法产品,令人扼腕。

对于前端拼接来说,要支持交互式全景类型如碗形、柱面等,同样也需要将一个全景播放模块嵌入平台软件,此种情况下,上述需求中除1、2、8外都仍然需要满足。

前端拼接技术的一个需要解决的问题是,以条带展开型全景作为全景帧类型做编码传输,如何节省带宽?将一个全景球展开为一个平面图像就如同将柚子皮拍扁在桌面上,总会像全球地图那样产生一个畸变,这是无法避免的。投影类型选择不好可能会导致相当大的畸变,比如球面两极地区一个像素被拉伸成一行像素。更好的投影类型应该是立方体展开

立方体展开

这一展开方式可以将畸变控制到很小的程度,但它一定程度上损失了条带全景图那种一览无余的直观性,需要特殊的全景播放器将它重新贴图到全景球或全景展开平面上才能够还原全景。当然也有其他采用更高级的数学方案设计的展开方式,这里不再提。

全景摄像机与全景直播

这里特意避免了提“VR全景”这一概念,因为严格来说VR全景和普通的全景摄像机并非同一概念,前者要求具有视觉深度感,后者只是个普通的2D曲面,沉浸感不强。但由于普通全景摄像机技术较前者简单,所以目前市面上大都为此类产品。

我的个人观点是,目前全景摄像机难以普及的一个关键是没有标准格式。并不像传统数码相机,全景输出格式杂乱无标准,全景视频播放器无法自动化决定采用何种投影类型播放,使得全景视频成了少数geek一族的玩具。但在真正的行业标准出现之前,让自己的产品对各种不同的输出格式都做好准备不失一个办法,而且不难。

这类产品中低端以前端拼的双鱼眼为主,高端以后端拼的多目摄像机为主,但迄今几乎没有很让人满意的产品出现。

双鱼眼方案的优势在于廉价且可以极小化拼缝。在所有可能的基于拼接算法的方案里,双鱼眼的拼缝是最小的。拼缝大小取决于多个摄像机投影中心的距离,摄像机的投影中心位置大致在sensor中心向后一个焦距远的地方,通常这是个很短的距离。理论上只有各个摄像机的投影中心重合于一点才能够产生出无缝的全景图,但这种情况下相机的体积需要压缩到极限,几乎不可达到,通常只能将尺寸压缩到极小以期更好的拼缝效果。除了双鱼眼方案,它是可以真的做到投影中心重合的。

所以,对于做基于拼接算法的全景摄像机的厂商,一个忠告是,将相机尺寸做小

全景直播机似乎有很长一段时间卡在很高的软件授权费和拼接服务器价格上,但这是比较奇怪的,因为这一技术并不困难——至少在安防领域,四年前就已经有公司做到了上千万像素的全景监控。像安防领域一样,最高性能且具有很好平台兼容性的方案就是OpenGL方案,现在的显卡处理几千万像素的全景拼接融合如同砍瓜切菜,顺便搞个硬编码做推流是不难的——我们自己的技术在这方面早已验证过。

全景直播机通常并不是多个摄像机拼一块儿这么简单粗暴,它需要解决两个基本问题,一是摄像机之间的帧同步,一是摄像机之间的成像参数同步。前者保证人通过拼缝时不会出现消失又出现这种诡异效果,后者使得全景画面亮度、色彩具有一致性,不出现尖锐的过渡。

但实际上,我们并不真的需要成像参数同步。理论上,多个摄像机各自自动曝光,可以实现HDR(高动态范围)全景,因而目前在硬件上做成像参数同步只是一个过渡方案,将来为了生成HDR效果全景,这一机制是必然要废弃的。

要有更好的拼接质量,可以选择CUDA或OpenCL,它比OpenGL提供更多控制力,使得开发者可以采用更复杂的图像处理算法。我们目前就在基于CUDA尝试HDR全景算法的开发。

VR全景

VR全景是万众期待众望所归,出于不可描述之原因,这一技术被视作新时代的宅男福利。但一定要冷静!因为我们真还有很多技术问题要解决。

实现3D效果,目前主要有基于传统拼接算法拼左右眼全景图(参见我们的文章《DIY 3D全景摄像机》)和光流算法(Google Jump/Facebook Surround360等)两种。

基于拼接算法基本是没有前途的(所以我们直接做成了DIY教程-_-!)。这一方案的死结在于,3D全景中深度感最强的近景,正好是拼缝最大的,而且你不能够通过缩小设备尺寸来解决,因为它至少应该有人的瞳距(~6.2cm)那么大,否则你戴上眼镜后,会发现自己缩小了——周边的一切都大了一遭。

光流算法是目前给出效果最好的,光流刻画了两个图像的像素是如何对应的,算法利用光流来插值计算没有被相机所采集的光线之颜色,从而产生出完全无缝的全景效果。但目前效率不高,关键是光流本身的计算是相当繁重的,而且算法对于每队图像还需要计算正反向两个光流,再考虑上光流在时间轴上的一致性,带来了非常大的计算开销。

实现VR视频采集,本质上是通过有限个相机采集几个点上的物理光线,然后用这些光线来猜测、插值出其他空间位置上任意光线。这在计算机视觉领域早已研究多年(想想黑客帝国里的子弹时间镜头是怎么来的),这个方向叫做”Image-Based Rendering”.

理想自然是通过采集有限个点上的光线就能够计算出一个邻域上的光场。这一定程度上做得到,而且有很好的工作,但付诸应用仍然有距离。

所以,仅就目前的情况来说,基于光流算法来做后期,做高质量近景VR视频是没问题的,但想要直播,还得等等。

统计算法的理论与实例之一 最大似然估计与EM算法

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

去年和朋友做了一些统计方面的东西,一些笔记发在这儿。统计我是外行,贻笑大方之处请指正。

\section{记号约定}
本文中用大写字母X,Y,Z...表示随机变量,小写字母x,y,z...表示随机变量的观测值,粗体表示向量/矩阵。\mathcal{X},\mathcal{Y},\mathcal{Z}...表示向量值随机变量或随机变量X,Y,Z...的观测样本。

\section{估计理论}
给随机变量一个待定参数的概率模型,通过一系列样本观测,萃取出概率模型最可能参数的理论。

\subsection{一些数学}
考虑随机变量X,Y的概率模型p(X,Y),我们经常需要考虑一系列观测样本\mathcal{X}=\{x_1,...x_N\},\mathcal{Y}=\{y_1,...y_N\}的分布

(1)   \begin{eqnarray*} p(\mathcal{X},\mathcal{Y})=\prod_i p(x_i,y_i) \end{eqnarray*}

观测样本之间的独立性意味着边缘分布p(\mathcal{X})(或p(\mathcal{Y})),条件分布p(\mathcal{X}\vert\mathcal{Y})(或p(\mathcal{Y}\vert\mathcal{X}))也满足类似的因子分解性质。严格说来,假设Y取值于\Delta,则\mathcal{Y}取值于\otimes\Delta^N,我们有

(2)   \begin{eqnarray*} p(\mathcal{X})&=&\sum_{\mathcal{Y}}p(\mathcal{X},\mathcal{Y})\\ &=&\sum_{\mathcal{Y}\in \otimes\Delta^N}\prod_i p(x_i,y_i)\\ &=&\prod_i\sum_{y_i\in\Delta} p(x_i,y_i)\\ &=&\prod_i p(x_i)\\ p(\mathcal{X}\vert\mathcal{Y})&=&\frac{p(\mathcal{X},\mathcal{Y})}{p(\mathcal{Y})}\\ &=&\prod_i \frac{p(x_i,y_i)}{p(y_i)}\\ &=&\prod_i p(x_i\vert y_i) \end{eqnarray*}

\subsection{最大似然估计}
假设随机变量X的概率模型为p(X\vert\theta),其中\theta为模型参数。我们有一系列独立观测样本\mathcal{X}=\{x_1,...x_N\},欲定出最可能的模型参数\theta_{\star}.

最大似然估计即假设观测结果是最大可能的结果,从而定出分布参数的估计方法。定义对数似然函数

(3)   \begin{eqnarray*} \log L(\theta;\mathcal{X})&=&\log p(\mathcal{X}\vert\theta)\\ &=&\log\prod_i p(x_i\vert\theta)\\ &=&\sum_{i}\log p(x_i\vert\theta) \end{eqnarray*}

则模型的最大似然估计为

(4)   \begin{eqnarray*} \theta_{\star}=\max_{\theta} \log L(\theta) \end{eqnarray*}

注意-\log p(x_i,\theta)是第i次观测所获取到的信息量,此目标函数一个含义是\textbf{调整参数使得观测到的信息量取极小}。

\subsection{概率分布的描述能力}
估计理论使得我们能够以估计分布的形式得到样本空间中的几何对象。概率分布的描述能力是广泛的,这里举几个例子。

\textbf{简单椭球区域}。可以简单的由Gauss分布p(\boldsymbol{X})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})刻画,\boldsymbol{\mu}给出了位置而\boldsymbol{\Sigma}则给出了椭球在各方向上的轴长。

\textbf{多个椭球区域}。由多个Gauss分布的叠加,即\textbf{混合Gauss分布}描述。概率密度分布\mathcal{M}(\{ \boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a,\boldsymbol{\pi}_a\})

(5)   \begin{eqnarray*} p(\boldsymbol{X})=\mathcal{M}(\{ \boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a,\boldsymbol{\pi}_a\}) \sim \sum_a \pi_a \mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a) \end{eqnarray*}

分布描述了样本空间里多个简单椭球区域,每个区域可以代表一个数据类。通常在应用中,我们把观测不到的隐含变量,类别标签Z加入到统计量集合中,于是概率密度函数成为

(6)   \begin{eqnarray*} p(\boldsymbol{X},Z)=\sum_a\delta_{Za}\pi_a\mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a) \end{eqnarray*}

\textbf{函数图形}。函数f(x)在样本空间里的图形可以由概率分布

(7)   \begin{eqnarray*} p(X,Y)=\mathcal{N}(f(X),\sigma^2) \end{eqnarray*}

刻画。给定任意X,Y都满足一个以f(X)为中心的Gauss分布,因而这里刻画了一个对Y=f(X)这一背后本质的数量关系所做的测量,而该测量是带有Gauss噪声的。

\subsection{最小二乘法}
考虑函数族f_{\boldsymbol{\alpha}}(x),我们通过一系列测量数据\mathcal{X}\equiv\{x_{1},...x_N\},\mathcal{Y}\equiv\{y_{1},...y_N\}去估计背后真实的函数关系f_*(x).
假设观测噪声为Gauss分布,于是观测的概率分布为

(8)   \begin{eqnarray*} p(X,Y)=\mathcal{N}(f_*(X),\sigma^2) \end{eqnarray*}

构造对数似然函数

(9)   \begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y})&=&\log p(\mathcal{X},\mathcal{Y}\vert\boldsymbol{\theta})\\ &=&\log\prod_i p(x_i,y_i\vert\boldsymbol{\theta})\\ &=&\sum_{i}\log \mathcal{N}(f_{\boldsymbol{\alpha}}(x_i),\sigma)\\ &=&\sum_{i} (-\log\sigma-\frac{1}{2}\log\pi -\frac{(y_i-f(x_i,\boldsymbol{\alpha}))^2}{2\sigma^2}) \end{eqnarray*}

于是,\max_{\theta} \log L(\theta;\mathcal{X},\mathcal{Y})给出

(10)   \begin{eqnarray*} &\max_{\boldsymbol{\alpha},\sigma}&\{-N\log\sigma-\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{2\sigma^2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha},\sigma}&\{N\log\sigma+\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{2\sigma^2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha}}&\{\frac{N}{2}\log\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{N}+\frac{N}{2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha}}&\{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2\} \end{eqnarray*}

其中\sigma=\sqrt{\sum_i(y_i-f_{*}(x_i))^2/N}.

于是我们得到了最小二乘法。

\subsection{EM算法}
对于一个联合分布,实际应用中并非每一个随机变量都可以被观测到。比如上述混合Gauss分布中,通常类别标签Z是无法观测到的。这种不完备观测情况下直接使用最大似然估计会出现计算量过大,实际不可解的问题。于是对于这类估计问题,我们使用\textbf{EM算法}来求解。

考虑混合Gauss分布

(11)   \begin{eqnarray*} p(X,Z\vert\boldsymbol{\theta})=\sum_a\delta_{Za}\pi_a\mathcal{N}(\mu_a,\sigma^2_a) \end{eqnarray*}

\mathcal{X}\equiv\{x_{1},...x_N\}为观测样本,\mathcal{Z}\equiv\{z_1,...z_N\}为观测样本对应的隐变量。对数似然函数为

(12)   \begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X})&=&\log\sum_{\mathcal{Z}} p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta})\\ &=&\log\sum_{\mathcal{Z}}\prod_i p(x_i,\mathcal{Z}\vert\boldsymbol{\theta}) \end{eqnarray*}

其中这一对数似然函数在计算上有诸多不便,首先加和号阻止了对数运算作用在p上,再者在z变量取值返回较广时计算代价也会很大。因此传统的最大似然估计在此处不适用。

在EM算法中,我们在每一个迭代步都构造一个新的函数Q^{(t)},用\max Q^{(t)}来近似\log L(\boldsymbol{\theta};\mathcal{X})在当前解\boldsymbol{\theta}^{(t)}附近邻域上的变化规律,迭代收敛到最终解

(13)   \begin{eqnarray*} Q^{(t)}(\boldsymbol{\theta})\equiv Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})&\equiv& E_{\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Z})]\\ \boldsymbol{\theta}^{(t+1)}&=&\mathrm{argmax}_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) \end{eqnarray*}

其中

(14)   \begin{eqnarray*} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)})= \frac{p(\mathcal{X},\mathcal{Z}\vert \boldsymbol{\theta}^{(t)})}{\sum_{\mathcal{Z}} p(\mathcal{X},\mathcal{Z}\vert \boldsymbol{\theta}^{(t)})} \end{eqnarray*}

该算法的正确性证明留到下一小节处理。

\subsection{EM算法的理论基础}
在目标函数较为复杂的情况下,极大/极小值的寻找通常采用迭代算法来处理。计算目标函数精确的形状是困难的,因而我们需要\textbf{在每一迭代步通过有限的计算,探知目标函数在当前迭代解之邻域上的变化规律,用一个简单的函数来刻画这一局域变化趋势,然后求解这一局部的极小/极大化问题,最后进入下一轮迭代}。这样,一个复杂目标函数的极小/极大问题就被分解成了一系列容易求解的简单问题,通过迭代不断靠近极优解。

具体说来,不同的算法会为目标函数的局域变化建立不同的模型,比如通过计算导数,用平面来建模该处的目标函数的变化;再如计算到二阶导数,用抛物面来近似目标函数,等等。

EM算法也是这样一种算法。\textbf{它的特殊之处在于,由于目标函数,即概然函数的特殊结构,我们可以找到比一次/二次近似更好的局域变化模型},从而获得比常规算法更好的性能。这一局域变化模型就是Q^{(t)}(\boldsymbol{\theta}).

首先我们来寻找Q^{(t)}(\boldsymbol{\theta})\log p(\mathcal{X}\vert \boldsymbol{\theta})之间的关系

(15)   \begin{eqnarray*} \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&\log p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta}) - \log p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta})\\ \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&\sum_{\mathcal{Z}} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)})\log p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta}) - \sum_{\mathcal{Z}} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \log p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta})\\ \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})+ H(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) \end{eqnarray*}

易证H(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})\geq 0,可见

(16)   \begin{eqnarray*} \log p(\mathcal{X}\vert \boldsymbol{\theta})&\geq&Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) \end{eqnarray*}

即\textbf{Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})构成对数似然函数的下界}。将15减去其在\boldsymbol{\theta}^{(t)}处的值,又有

(17)   \begin{eqnarray*} \Delta \log p(\mathcal{X}\vert \boldsymbol{\theta}) - \Delta Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) &=& H(\boldsymbol{\theta}\vert \boldsymbol{\theta}^{(t)}) - H(\boldsymbol{\theta}^{(t)}\vert \boldsymbol{\theta}^{(t)})\\ &=&D(p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \mid p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}))\\ &\geq& 0 \end{eqnarray*}

即\textbf{每一迭代步中必有\Delta \log L \geq \Delta Q}.

可见Q^{(t)}(\boldsymbol{\theta})可以反映\log p(\mathcal{X}\vert \boldsymbol{\theta})\boldsymbol{\theta}^{(t)}处的变化规律。更重要的是,Q^{(t)}(\boldsymbol{\theta})构成了对数概然函数的一个下界,这样以来我们就可以放心的求解\max Q^{(t)}(\boldsymbol{\theta})这一问题来得到一个更好的迭代解,而不用像LM算法中那样担心二阶模型只在一个局域的范围内构成目标函数的有效近似,需要用置信域方法避免大步长的优化。

算法的收敛性由这一点来保证。由于D(p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \mid p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}))\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}附近的一阶展开项为零,于是17的一阶展开给出

(18)   \begin{eqnarray*} \frac{\mathrm{d}}{\mathrm{d} \boldsymbol{\theta}} \log p(\mathcal{X}\vert \boldsymbol{\theta})\vert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}} &=& \frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\theta}} Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})\vert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}} \end{eqnarray*}

算法迭代到Q不再有明显变化(\boldsymbol{\theta}^{(t+1)}\approx\boldsymbol{\theta}^{(t)})时,上式右端为零,式18意味着\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}同样也是对数似然函数的极值点。于是Q的收敛意味着对数似然函数的收敛,算法必然收敛于一个极优解。

Notes on Conjugate Gradient Method – 2

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

打算再将这个主页做一做,先更一波旧笔记。还是那句话,这不是教程,不适合初学。本文中记号大致与前一篇Notes on Conjugate Gradient Method一致。

\mat{A}对称正定时,共轭梯度法是一个很好的选择。

此时的极小化问题其实是一个\mat{x}\mat{x}^{*}间的距离极小化问题

(1)   \begin{eqnarray*} &&\min_{\mat{x}}\frac{1}{2}\mat{x}^{T}\mat{A}\mat{x}-\mat{b}^{T}\mat{x}\\ \Leftrightarrow&&\min_{\mat{x}}\lVert \mat{x}-\mat{x}^{*} \rVert_{\mat{A}} \end{eqnarray*}

\mat{A}为度量,在一个\mat{A}正交标架下,上述距离极小化问题成为一个简单的欧式距离极小化问题,只需在依次各\mat{A}正交方向上极小化距离即可。
这需要解决以下问题:
\begin{enumerate}
\item 一个搜索方向上如何极小化\lVert \mat{x}-\mat{x}^{*} \rVert_{\mat{A}}
\item 如何高效的构造出与之前搜索过的子空间\mat{A}正交的新搜索方向;
\end{enumerate}

第一个问题。由简单的几何直观,距离的极小化等价于

(2)   \begin{eqnarray*} &&\langle\mat{e}^{(n)},D_{n}\rangle_{\mat{A}}=0\\ \Leftrightarrow&&\langle\mat{r}^{(n)},D_{n}\rangle=0 \end{eqnarray*}

其中D_{n}=\mathrm{span}\{ \mat{d}^{(0)},\dots,\mat{d}^{(n-1)} \}表示\mat{x}=\mat{x}^{(n)}时已经搜索过的子空间。

本式D_{n}换成\mat{d}^{(n)}同样成立。所以,只需搜索到\mat{r}^{(n)}正交于搜索方向\mat{d}^{(n)}时为止即可达到极小化距离的目的。

第二个问题。我们令D_{n}=\mathcal{K}_{m}(\mat{A},\mat{r}^{(0)}),其中\mat{r}^{(0)}为初始残差,于是有D_{n}=\mathrm{span}\{ \mat{r}^{(0)},\mat{A}D_{n-1} \}.这样式2就意味着

(3)   \begin{eqnarray*} \langle\mat{r}^{(n)},D_{n-1}\rangle_{\mat{A}}=0 \end{eqnarray*}

这带来一个极大的便利,即算法可以通过\mat{r}^{(n)}快速构造出与D_n子空间\mat{A}正交的\mat{d}^{(n)},为此只需要将\mat{r}^{(n)}\mat{d}^{(n-1)}方向的分量(\mat{A}内积意义上)扣除即可。

这样我们可以得到共轭梯度法的算法流程:每次搜索都搜索到残差方向与搜索方向正交为止,新的搜索方向由当前残差用\mat{A}内积与\{ \mat{d}^{(0)},\dots,\mat{d}^{(n-1)} \}做正交化得到,实际上只需要扣除\mat{d}^{(n-1)}方向的分量即可。
算法描述为

(4)   \begin{eqnarray*} \mat{x}^{(n+1)}&=&\mat{x}^{(n)}+\alpha^{(n)}\mat{d}^{(n)}\\ \mat{d}^{(n)}&=&\mat{r}^{(n)}-\beta_{n,n-1}\mat{d}^{(n-1)} \end{eqnarray*}

其中\alpha\mat{r}^{(n+1)}\mat{d}^{(n)}的正交性给出

(5)   \begin{eqnarray*} &&\mat{r}^{(n+1)^T}\mat{d}^{(n)}=0\\ \Leftrightarrow&&(\mat{e}^{(n)}+\alpha_{n}\mat{d}^{(n)})^T\mat{A}\mat{d}^{(n)}=0\\ \Leftrightarrow&&\alpha_{n}=-\frac{\langle \mat{d}^{(n)},\mat{e}^{(n)} \rangle_{\mat{A}}}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}=\frac{\langle \mat{d}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}\\ &&=\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}} \end{eqnarray*}

\beta则这样计算出

(6)   \begin{eqnarray*} &&\mat{r}^{(n)}=-\mat{A}\mat{e}^{(n)}=\mat{r}^{(n-1)}-\alpha_{n-1}\mat{A}\mat{d}^{(n-1)}\\ \Leftrightarrow&&\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle = -\alpha_{n-1} \langle \mat{r}^{(n)},\mat{d}^{(n-1)} \rangle_{\mat{A}}\\ \Rightarrow&&\beta_{n,n-1}=-\frac{\langle \mat{r}^{(n)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}{\langle \mat{d}^{(n-1)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}\\ &&=\frac{1}{\alpha^{(n-1)}}\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n-1)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}\\ &&=\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{r}^{(n-1)},\mat{r}^{(n-1)} \rangle} \end{eqnarray*}

共轭梯度法在加入precondition时会遇到一个问题,即如何保证\mat{M}^{-1}\mat{A}仍是对称矩阵。
对于对称正定矩阵\mat{M},存在一个分解\mat{M}=\mat{E}\mat{E}^T.容易证明\mat{M}^{-1}\mat{A}\mat{E}^{-1}\mat{A}\mat{E}^{-T}具有相同的本征值,因而具有同样的precondition效果。于是可以将问题转化为

(7)   \begin{eqnarray*} \mat{E}^{-1}\mat{A}\mat{E}^{-T}\hat{\mat{x}}=\mat{E}^{-1}\mat{b},\hat{\mat{x}}=\mat{E}^{T}\mat{x} \end{eqnarray*}

再用CG求解。\mat{E}^{-1}\mat{A}\mat{E}^{-T}的几何意义是将问题转入一个新的坐标系下求解,这个坐标系下二次型的超椭球被变换的更接近超球了一些。这样要求计算出\mat{E},但实际上我们可以不用计算\mat{E}而直接用\mat{M}完成计算。为此只需要将原来的迭代公式做坐标变换,写出新坐标系下的迭代公式即可。要注意迭代公式在该变换下并不完全保持形式一致,具体的对应关系如下所示

Rendered by QuickLaTeX.com
有稍许微分几何或相对论背景的同学不难看出这其实是一组在坐标系变换下具有不同变换性质的量之变换式,包含了协/逆变矢量,标量的变换性质,以及如何将当前坐标系下定义,非形式不变的量\langle \mat{r},\mat{r}^{\prime} \rangle提升为任意坐标系下均保持形式不变的量\langle \mat{r},\mat{M}^{-1}\mat{r}^{\prime} \rangle的操作.

替换后可将\mat{E}全部消除掉,剩下仅用\mat{M}表达的迭代公式

(8)   \begin{eqnarray*} \mat{d}^{(0)}&=&\mat{M}^{-1}\mat{r}^{(0)}\\ \mat{d}^{(n)}&=&\mat{M}^{-1}\mat{r}^{(n)}-\beta_{n,n-1}\mat{d}^{(n-1)}\\ \mat{x}^{(n+1)}&=&\mat{x}^{(n)}+\alpha^{(n)}\mat{d}^{(n)}\\ \alpha_{n}&=&\frac{\langle \mat{r}^{(n)},\mat{M}^{-1}\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}\\ \beta_{n,n-1}&=&\frac{\langle \mat{r}^{(n)},\mat{M}^{-1}\mat{r}^{(n)} \rangle}{\langle \mat{r}^{(n-1)},\mat{M}^{-1}\mat{r}^{(n-1)} \rangle} \end{eqnarray*}

DIY 3D全景摄像机

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

目前市面上的VR全景多是二维,没有深度感,若真想拥有身临其境般的体验,深度信息必不可少。诺基亚的OZO,Google的Jump,Facebook开源的Surround 360,都是为3D全景而设计。OZO设备8个鱼眼售价高达三十多万,Surround 360搭载的是Point Grey的相机,硬件成本二十多万,Jump也要搭载十几个GoPro,硬件成本少说也要几万,普通玩家真心想玩也要思量的。如果不追求那么高大上,其实自己就可以DIY出一台3D全景相机。

组件:

千兆交换机+网络摄像头模组+180°鱼眼镜头+线材

结构:

要实现深度感,结构是关键。

相机可以理解为光线采集设备,采集到的光线与成像平面的交点即像点。

通常的二维全景要求采集到的所有光线汇于中心点,即视点,以视点为中心的球面或圆柱面为成像面,所有光线交于成像面形成全景图,如下图(a)。二维全景相机要求所有相机共中心摆位,即所有相机的光轴相交于视点。

3D全景即左右眼各对应一个全景图。

两只眼睛分别对应两个不同的视点位置,当转头360度时,两只眼睛转过的轨迹即一个以瞳距为直径的圆,称之为Viewing Circle,3D全景要求采集到的所有光线相切于此圆。左右眼采集到的光线分别与成像面相交形成左眼全景图和右眼全景图, 如下图(b)(c)。

两只眼睛所在的视点位置投影出的图像称为一个立体对,如图(b)左眼光线1和图(c)右眼光线1即可看成一个立体对,同理左眼中的光线2,3,4,5等与右眼中的2,3,4,5等分别构成不同的立体对。

该图引自文献:Stereo Panorama with a Single Camera.

至于相机摆位一般有两种方案,如下图:

该图引自文献:Jump: Virtual Reality Video.

切向摆位如上图(a):每个相机的光轴相切于Viewing Circle,此种方案一半的相机用于左眼全景(图(a)中绿色相机),而另一半的相机用于右眼全景(图(a)中红色相机)。

径向摆位如上图(b): 每个相机的光轴沿Viewing Circle半径方向,此种方案不区分左右眼相机,每个相机都对两眼的全景图有贡献,因此对每个相机水平视场角有更高要求:R越小,要求每个相机的水平视场角越大。

一般R设计比较大时采用径向摆位,R较小时采用切向摆位。

切向摆位最简单的结构设计即正多边形,每条边上放置两个camera,其sensor中心的距离设为瞳距。如果镜头的视场角足够大,可以设计一个正三角形,用六个camera来实现3D全景。本文介绍的是正四边形八个camera的方案,用Solidworks设计一个简单的支架,预留出上camera的安装孔位。

结构设计及3D打印: 

组装:

效果:

原始视频截图

上下3D格式

VR眼镜观看3D效果

目前能实现3D全景的技术无外乎几种:

  • 拼接方案。左右眼的视野分别做拼接融合以达到3D全景的效果,这种方案最简单,可以实时化,其缺点是拼缝难消除。
  • 光流方案。Google Jump以及Facebook开源的算法即此方案,能很好的消除拼缝,但实时化比较困难,适合做后期处理。
  • 光场重建。用有限个相机重建光场,给出真正的3D效果,使用户拥有更多活动自由度,这是终极的VR视频采集方案,但即便在理论上也有很多困难之处。

本文中所用的实时3D全景拼接软件是奇点视觉实时全景拼接方案,演示中为2700万像素实时拼接融合效果。我们目前可以轻易的实现超高分辨率实时拼接和直播,但仍有以下几个问题:

  1. 硬件使用了安防用网络摄像头模组,不具备同步曝光功能,因此全景图明暗不均较严重;
  2. 近景拼缝明显。3D全景要求设备直径至少等于人眼瞳距,但过大的直径容易导致更严重的拼缝。这是基于拼接方案做3D全景的一个本质困难。

但至少到目前,拼接方案仍然是唯一能做到低开销高分辨率直播的。我们正在研发新一代的实时光流/光场重建算法,希望能解决高质量近景VR直播的问题。

奇点视觉是一个致力于计算机视觉技术的研发和产品化的团队,专注于算法,为了发挥我们的优势,没有去做产品,而是给有能力做产品的公司提供技术解决方案。过去两年,我们专注于安防全景,毫不夸张的说,安防全景技术我们做到了世界顶级水平,已经由客户厂商实现产品化并销往海内外。目前我们致力于把全景技术迁移到VR应用中,并作技术升级,希望能够为更高质量的VR内容贡献一份力量,敬请关注奇点视觉。

全景视频相关需求征集

首先说一下我们目前的情况。到今天,我们的全景技术已经非常完善的产品化,性能、质量和稳定性都趋于收敛,可以不夸张的说,这一技术目前处于世界上领先的地位。

我们有目前最高的实时拼接融合性能,在普通笔记本上可以做到数千万像素的实时拼接。在这一实时拼接模块的驱动下,实现超高清全景直播只需要一台普通中端笔记本。基于这一模块,我们在安防领域实现了4000w像素实时全景监控,且在笔记本上,子码流模式下可以同时进行多达30个全景实时浏览。该模块其余特征包括:

  1. 子/主码流动态切换机制。小窗口时用子码流,放大观看用主码流,后台实时切换,非常适合安防监控领域。
  2. 多种全景类型。包括但不限于交互式球面/半球/柱面,球面180/360展开,柱面180/360展开,小世界,双目上下/左右格式的3D全景等等。
  3. 可以将全景图实时播放,也可以实时回传给图像分析算法。
  4. 可以将全景像素坐标、输入视频像素坐标、物理空间实际方向互相转换。在安防领域可以以此实现全景/球机联动,全景/全景联动实现局域放大等;也可以籍此实现一些增强现实效果。

我们有一套强大的量产标定算法,支持对任意多路,任意多不同类型的摄像机进行最优拼接参数标定,支持无显著重合区情况下的标定,拼缝达到理论极小值,固定步骤可以得到稳定标定结果,返工率极低

以上全部技术完全自主知识产权。不依赖于PTGUI/OpenCV等第三方实现。也正因为如此,我们对技术中全部细节都具有掌控力,不会因为一个问题出在第三方代码中而束手无策。

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

然后说一下此文的真正目的,即收集各行业对于全景视频的需求

VR这个风口吹起来了漫天的塑料袋,但没有看到太多扎实漂亮的产品。我们希望做一款这样的产品。所以我们希望收集您对于全景视频最迫切的渴求,以此决定接下来的方向。您只需简略说明在您这里全景的用途,以及对它的一些特殊要求即可,可以在本站留言,也可发邮件:

planckscale1729@163.com

QQ是397692433, 但不太经常上,可能没法及时相应您:-P

升维成功,有3D全景了

最近出门,在深圳的旅馆里埋头干了两天,算法终于可以支持3D全景了,昨晚十一点多用CardBoard看到自己第一个3D全景视频,很幸福。

得益于初期花费大量精力建立起来的良好架构,加入3D全景的过程还算轻松,另外由于近景标定技术的支持,我们可以对大半径的3D全景相机做高精度标定,达到理论最优效果。早期投入的时间和精力还算没有白费。

3D全景技术目前大致有两类,一类基于传统的拼接,双目各自拼一个全景球;一类基于双目的三维重建。我们目前是第一类技术。

第一类技术门槛较低,是目前市面上能见到的主要方式。这类技术的缺陷主要是拼缝难以克服,另外理论上无法实现真正的3D全景,在抬头低头时没有深度感。具体的体验能达到怎么样的程度,也是我们接下来要重点测试的。

采用第二类技术的典型是google jump, 目前据说jump仍然要消耗巨量的计算力来完成这种重建,直播更是不可能。它的好处在于可以实现真正的3D效果,而且无缝。更快的算法应该并非不可能,而且很可能已经有公司在致力于这件事情。

话说最近VR Porn很火,谁看谁知道。

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

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

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

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

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

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

QQ: 397692433

手机: 15165701250

全景摄像机量产标定算法

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

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

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

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

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

QQ: 397692433

手机: 15165701250

高清全景视频拼接技术 简介

目前为止我们的全景技术有三种:

  1. 高质量融合全景。基于普通摄像机,用CUDA显卡进行实时拼接和高质量融合的拼接技术。基于该种技术,可以用相对较差的摄像机输出非常好的全景视频,对相机的相对摆位要求最为宽松。但这一技术对计算机性能有一定要求,需要有CUDA显卡(Nvidia近几年的显卡都支持CUDA技术)的支持。该技术详情可见《多路视频实时全景拼接算法》一文的介绍。
  2. 高效融合全景。这种技术采用了相对简单的融合算法,可以在极低的开销下产生出高分辨率全景输出,在各种硬件平台(中低端台式机、笔记本,甚至移动平台)上都能低开销流畅运行。由于融合算法简单,因此要产生高质量全景输出,全景摄像机就需要做到曝光、增益等参数的同步。
  3. 鱼眼全景。这类技术最简单,基于鱼眼镜头,将鱼眼画面矫正为更易浏览的平面展开,或投影到全景球面上做沉浸式浏览。

第一种技术之前已经有文章着重介绍,接下来我们将会发布关于后两种技术的一系列Demo视频以及性能数据。

基于高效融合技术的8路960p视频实时全景浏览

不是意外的意外,全景视频技术在最近几个月爆发式的进入公众视野。抛掉概念炒作的浮云不说,这项技术至少在安防领域将会有不小的市场。作为还算是稀有动物的全景视频核心技术提供商之一,我们不应该沉默下去了。

我们今后原则上不再销售SDK,而主要以技术合作为主。欢迎有具备硬件研发能力,对做高质量、高水准产品有兴趣的同道来电洽谈合作。

QQ: 397692433
Phone: 15165701250
Email: planckscale1729@163.com