作者归档：steven

update.

或许熟悉我们的朋友已经知道，我们大概有两年没有在全景技术方面继续展开新的工作了。这里的主要原因是，我们跳进了另一个坑，算法交易，并且作为一个袖珍团队，很快就被吸干了精力。另一方面，像全景这样的技术，通过技术合作或者技术服务来盈利并非易事，你很难确保知识产权不被侵犯（企业跟你边“合作”边“学习”，事了另起炉灶留你干瞪眼儿这种事儿清华计算机系的老师们都没少遇到），另外整个产品化周期过长，自身团队之外的种种不可控因素过多，都成为要考虑的问题。最初选择这种形式来运营只有一个原因，就是自由和独立。保持团队有限精力不被过多侵蚀，保持做“闲事儿”的权利。如今跳了另一个坑，也无非是这样一个动机。

由于两年间精力都在构筑我们的算法交易平台，我并没有太多时间来考虑全景技术的未来。直到后来看到海思推出了芯片解决方案，我心里一沉，感觉14年时就曾经在担心的事情，终于还是来了。去年我开始考虑将一部分未曾做落地应用，单纯玩儿出来的算法做开源，但并没有足够多时间来整理这些工作，迄今仍未实现。

我们经常嘲笑自己，埋头在家做了好多东西，很多都是工业级强度的东西，都扔在硬盘里睡大觉，就又跑去搞别的新奇玩意儿。这是一个纯技术团队的很大缺陷，而我们直到今天也没（lan）有（de）克服。

所以今天对于我们全景技术的未来做一些说明。

首先，之前的合作我们一直在持续，产品也一直都稳定量产着，有兴趣的客户可以联系深圳华途数字了解产品具体情况。这些高清晰度安防全景产品应用了我们的量产标定技术以及实时拼接融合技术，由我们提供二次开发SDK,算法质量顶级（量产三年来大概只fix了1次较严重bug[在某些型号的机器上，某个全景投影类型和别的类型混淆，已经完全通过软件升级解决]，加上我不记得的小bug，总量应该小于5个），您可以信赖。

然后，我们会尽可能花空闲时间整理过去的未公开技术并完成开源。至于工业级代码是否开源，取决于今后的具体情况。

最后，目前已经产品化的量产标定技术、实时拼接融合技术，如果您是安防或VR业内top5的靠谱公司，欢迎来聊聊技术收购（客官我们只卖术不卖身吆）。如果您是海思这样的高富帅芯片厂商，请直接来谈朋友（求你了）！如果你们都不来，我就择吉日都开源掉，我们这样的geek还是适合做开源。

至于这个网站，以后会自我放飞成为一个geek的后花园。

统计算法的理论与实例之二 EM算法实例： 形状约束的非参数混合高斯模型

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

最近比较无聊，干脆把这个系列补完。本文的内容是我表妹之前在芝大读书时的一个小题目，作为EM算法的一个比较复杂的实例来演示是不错的。由于我没找到引文，就不引用了。我不是做统计的，有错误欢迎留言。

\subsubsection{解析推导}
这是一个更为复杂的估计问题实例。要估计的分布是样本空间里一组函数图形对应的分布。具体说来，我们有随机变量\{\boldsymbol{X},Y,Z\},其中\boldsymbol{X}\in \mathbf{R}^p,Y\in \mathbf{R},Z \in \{1,...K\},其中Z是混合高斯模型中的类别标签变量,对于不同的Z值,Y和\boldsymbol{X}存在一个不同的潜在函数关系m_Z(x),其中m_Z(x)为凸函数。概率模型表述为

(1)   \begin{eqnarray*} p(\boldsymbol{X},Y,Z)&=&\sum_j\delta_{Zj}\pi_j\mathcal{N}(m_j(\boldsymbol{X}),\sigma^2_j) \end{eqnarray*}

现在假定Z不可观测，我们从一系列观测值\{\boldsymbol{x}_i,y_i\}中估计出模型的参数\boldsymbol{\theta}\equiv \{\boldsymbol{\pi},\boldsymbol{m(x)},\boldsymbol{\Sigma}=\{\sigma_1,...\sigma_{K}\}\}的真实值\boldsymbol{\theta}_*。

接下来我们用EM算法来处理这一问题。设\mathcal{X}\equiv\{\boldsymbol{x}_{1},...\boldsymbol{x}_N\},\mathcal{Y}\equiv \{y_1,...y_N\}为观测样本,\mathcal{Z}\equiv\{z_1,...z_N\}为观测样本对应的隐变量。我们首先计算p(\mathcal{Z}\vert \mathcal{X},\mathcal{Y},\boldsymbol{\theta}^{(t)})

(2)   \begin{eqnarray*} \gamma_{ij}^{(t)}&\equiv&p(z_i=j\vert \boldsymbol{x}_i, y_i, \boldsymbol{\theta}^{(t)})\\ &=&\frac{p(\boldsymbol{x}_i,y_i,z_i=j \vert \boldsymbol{\theta}^{(t)})}{\sum_k p(\boldsymbol{x}_i,y_i,z_i=k \vert \boldsymbol{\theta}^{(t)})}\\ &=&\frac{\pi_j \mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}{\sum_k \pi_k \mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}} \end{eqnarray*}

类似于边缘分布和条件分布的因子分解性质，可以证明对任意满足q(\mathcal{Z})=\prod_i q_i(z_i)的概率分布q(\mathcal{Z}),有

(3)   \begin{eqnarray*} E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]=\sum_i E_{q_i(z_i)}[\log L(\boldsymbol{x}_i,z_i)] \end{eqnarray*}

这意味着Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})可以分解成在每个样本上各自均值之总和的形式

(4)   \begin{eqnarray*} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) &\equiv& E_{\mathcal{Z}\vert \mathcal{X},\mathcal{Y},\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y},\mathcal{Z})]\\ &=& \sum_i E_{z_i\vert \boldsymbol{x}_i,y_i,\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i)] \end{eqnarray*}

我们把该性质的证明留到后面。

于是

(5)   \begin{eqnarray*} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) &=&\sum_{ij} p(z_i=j\vert \boldsymbol{x}_i,y_i,\boldsymbol{\theta}^{(t)}) \log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i=j)\\ &=&\sum_{ij} \gamma_{ij}^{(t)} \log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i=j)\\ &=&\sum_{ij} \gamma_{ij}^{(t)} [\log \pi_j - \log\sigma_j - \frac{1}{2}\log{2\pi} - \frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2}] \end{eqnarray*}

需要注意的是\boldsymbol{m}(\boldsymbol{x})=\{\boldsymbol{m}_1(\boldsymbol{x}),...\boldsymbol{m}_{K}(\boldsymbol{x})\}包含了无穷个参数需要优化，为了约化这一问题，我们将其离散化，取m_{ji}=m(\boldsymbol{x}_i)为优化参数，仍然简记为\boldsymbol{m}.这样，\boldsymbol{m}(\boldsymbol{x})被约化到用K\times N个参数来表示。

接下来我们需要求解的极大化问题为

(6)   \begin{eqnarray*} \max_{\boldsymbol{\theta},\boldsymbol{\beta}} &&Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})\\ \rm{subject\ to} &&\sum_j \pi_j=1\\ &&m_{ji}-m_{jk} \geq \boldsymbol{\beta}_{jk}^{T} (\boldsymbol{x}_i-\boldsymbol{x}_k)\qquad \forall i,j,k \end{eqnarray*}

这里我们为m(x)加入了凸函数的形状约束。约束的几何意义很明确，即\textbf{函数总是位于每一点的切线之上方}。

由于Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})中\boldsymbol{\pi}的极大化和其他参数相互独立，因而可将问题解耦。问题简化为两个子问题，其中\boldsymbol{\pi}参数对应的子问题为

(7)   \begin{equation*} \left\{ \begin{aligned} \boldsymbol{\pi}^{(t+1)}&=argmax_{\boldsymbol{\pi}} \sum_{ij}\gamma_{ij}^{(t)}\log\pi_j\\ \sum_j \pi_j&=1 \end{aligned} \right. \end{equation*}

可以求出解析解，即

(8)   \begin{eqnarray*} \pi_j^{(t+1)}&=&\frac{\sum_{i}\gamma_{ij}^{(t)}}{\sum_{ij} \gamma_{ij}^{(t)}}=\frac{\sum_i \gamma_{ij}^{(t)}}{N} \end{eqnarray*}

\boldsymbol{m}参数对应的子问题

(9)   \begin{equation*} \left\{ \begin{aligned} \boldsymbol{m}^{(t+1)}&=argmax_{\boldsymbol{m},\boldsymbol{\beta}} \sum_{ij}\gamma_{ij}^{(t)}[- \log\sigma_j - \frac{(y_i-m_{ji})^2}{2\sigma_j^2}]\\ &=\sum_j argmax_{\boldsymbol{m},\boldsymbol{\beta}}[-N\log\sigma_j-\frac{\sum_i \gamma_{ij}^{(t)}(y_i-m_{ji})^2}{2\sigma_j^2}]\\ &m_{jk}-m_{ji} \geq \boldsymbol{\beta}_{ji}^{T} (\boldsymbol{x}_k-\boldsymbol{x}_i)\qquad \forall i,j,k \end{aligned} \right. \end{equation*}

是个二次规划问题。由于\sigma_j,m_{ji}只出现在具有相同j值的项中，这意味着按不同的j值把项分组，每一组都是独立变化的，因而可以作为一个独立的子问题来最大化。可以发现，\textbf{这里的Gauss混合模型中每个component都对应于一个带约束的最小二乘问题}。这一性质可以用来降低计算复杂度，提高算法效率。

在一些情况下，二次规划中的约束可以大幅度的约减。第一种情况是当X为单变量时。不失一般性的，我们假设观测样本已经按照X的大小做升序排列，此时可以用二阶导数非负来作为凸性约束。于是约束可约化为

(10)   \begin{eqnarray*} m_{jk}-m_{ji} &\geq& \beta_{ji} (x_{k}-x_i)\qquad \forall i,j,k\in\{i-1,i+1\}\\ \beta_{j,i+1}-\beta_{j,i}&\geq&0 \qquad\forall i\in\{1...N-1\},j \end{eqnarray*}

这里\beta_{ji}代表了第j个component在第x_i处导数的估计值。需要注意的是，为了得到一点上导数（切线）的估计，我们需要两个约束，约束在左右两侧的函数值都要位于切线上方。

另一种情况则是当X为多维的向量，且有m(\boldsymbol{X})=\sum_\alpha m^\alpha(X_\alpha)时。此时每个分量m^\alpha(X_\alpha)的凸性意味着m(\boldsymbol{X})也是凸的，因此约束可写为

(11)   \begin{eqnarray*} m^{\alpha}_{jk}-m^{\alpha}_{ji} &\geq& \beta^{\alpha}_{ji} (x_{k}-x_i)\qquad \forall \alpha,i,j,k\in\{i-1,i+1\}\\ \beta^{\alpha}_{j,i+1}-\beta^{\alpha}_{j,i}&\geq&0 \qquad\forall \alpha,i\in\{1...N-1\},j \end{eqnarray*}

即对向量函数的每个分量，都满足\eqref{mConstraint1D}所示的约束。

最后我们来证明E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]的分解性质，其中q(\mathcal{Z})为任意满足q(\mathcal{Z})=\prod_i q_i(z_i)的概率分布

(12)   \begin{eqnarray*} &&E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]\\ &=&\sum_{\mathcal{Z}} q(\mathcal{Z})\log p(\mathcal{X},\mathcal{Z})\\ &=&\sum_i \sum_{\mathcal{Z}} q(\mathcal{Z}) \log p(\boldsymbol{x}_i,z_i) \end{eqnarray*}

可见对每个sample，都需要在分布q(\mathcal{Z})上求\log p(\boldsymbol{x}_i,z_i)的均值。但由于\mathcal{Z}中包含大量无关的z_k\neq z_i,它们都应该被加和为1,从而只剩下q_i(z_i):

(13)   \begin{eqnarray*} &=&\sum_i \sum_{\{z_1,...,z_N\}} \prod_j q_j(z_j) \log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i \sum_{z_i}\{\underbrace{\sum_{\{z_1,.\slashed z_i.,z_N\}} \prod_{j\neq i} q_j(z_j)}_{=1}\} q_i(z_i) \log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i \sum_{z_i} q_i(z_i)\log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i E_{q_i(z_i)}[\log L(\boldsymbol{x}_i,z_i)] \end{eqnarray*}

\subsubsection{算法实现}
算法实现主要需要注意的是防止数值下溢带来的除零问题。技术上，避免概率连乘积，避免指数，让计算的中间结果尽量在对数域中表示是一些基本技巧。

首先，在计算\gamma_{ij}时，可能会出现某个样本点远离所有component的情况，此时分子分母都容易下溢为零，而出现\frac{0}{0}的情况。我们做如下处理

(14)   \begin{eqnarray*} \gamma_{ij}^{(t)}&=&\frac{\pi_j \mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}{\sum_k \pi_k \mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}}\\ &=&\frac{1}{\sum_k \frac{\pi_k}{\pi_j} \frac{\mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}}{\mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}}\\ &=&\frac{1}{\sum_k \frac{\pi_k}{\pi_j} \frac{\sigma_j}{\sigma_k} \exp(-\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}+\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})}\\ \end{eqnarray*}

算法中需要计算对数似然函数的上升率，来判断算法的终止条件是否满足。对数似然函数中有大量的概率连乘积，我们作如下处理来消除之

(15)   \begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y}) &=&\log p(\mathcal{X},\mathcal{Y}\vert\boldsymbol{\theta})\\ &=&\log\prod_i p(\boldsymbol{x}_i,y_i\vert\boldsymbol{\theta})\\ &=&\sum_i\log\sum_{z_i} p(\boldsymbol{x}_i,y_i,z_i\vert\boldsymbol{\theta})\\ &=&\sum_i\log \sum_j \frac{\pi_j}{\sigma_j \sqrt{2\pi}} \exp (-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\\ &=&\sum_i\log \frac{\pi_k}{\sigma_k \sqrt{2\pi}} \exp (-\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2})\\ && \sum_j \frac{\pi_j}{\pi_k}\frac{\sigma_k}{\sigma_j} \exp (\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\\ &=&\sum_i\{\log \frac{\pi_k}{\sigma_k \sqrt{2\pi}} -\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}\\ && +\log\sum_j \frac{\pi_j}{\pi_k}\frac{\sigma_k}{\sigma_j} \exp (\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\} \end{eqnarray*}

对于二次规划子问题，我们需要将其化为标准形式，从而可以利用现成的二次规划库来求解。

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

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(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*}

全景视频相关需求征集

发表回复

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

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

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

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

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

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

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

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

planckscale1729@163.com

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

开卖了

发表回复
Shenzhen Civic Center

Shenzhen Civic Center

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

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

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

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

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

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

5 6 13 14

升维成功，有3D全景了

2条回复

最近出门，在深圳的旅馆里埋头干了两天，算法终于可以支持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