Author Archives: steven

CUDA, 软件抽象的幻影背后 之二

先更新到这儿,稍后再回来抛光查错。CUDA比较杂,我一写起来容易满嘴跑火车弄出错误,欢迎拍砖。

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

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

上一篇里说到,有两点对CUDA的计算能力影响甚大:数据并行,以及用多线程掩盖延迟。接下来我们要深入到其硬件实现,看一看这些机制是如何运作的。

通常人们经常说某GPU有几百甚至数千的CUDA核心,这很容易让人联想到多核CPU。不过事实上两种“核心”是不一样的概念,GPU的CUDA核心只相当于处理器中的执行单元,负责执行指令进行运算,并不包含控制单元。可以类比到CPU核心的是流多处理器(Streaming Multiprocessor,简写为SM. Kepler中叫做SMX,Maxwell中叫做SMM),通常一个GPU中有数个SM,而每个SM中包含几十或者上百个CUDA核心,以及数个warp scheduler(相当于控制单元)。如下图GM204中有16个SM,每个SM中有128个CUDA核心,4个warp scheduler。

GeForce_GTX_980_SM_Diagram-545x1024

图 1.  GM204的SM结构图

每个SM中有大量的寄存器资源,在GM204的例子中,有总共64k 32-bit寄存器,可以养活成千上万的线程。SM中另外一个重要资源是Shared Memory,没错,它正是软件抽象中Shared Memory的对应物。在GM204中,每个SM有96KB的Shared Memory.

到这里,SM在软件抽象里的对应也呼之欲出了,没错,正是Block。我们不妨先摆出这个对应:
Block <-> SM
Thread执行 <-> CUDA Cores
Thread数据 <-> Register/Local Memory

同一Grid下的不同Block会被分发到不同的SM上执行。SM上可能同时存在多个Block被执行,它们不一定来自同一个kernel函数。每个Thread中的局域变量被映射到SM的寄存器上,而Thread的执行则由CUDA核心来完成。

SM上可以同时存在多少个Block?这由硬件资源的消耗决定:每个SM会占用一定数量的寄存器和Shared Memory,因此SM上同时存活的Block数目不应当超过这些硬件资源的限制。由于SM上可以同时有来自不同kernel的Block存在,因此有时候即便SM上剩余资源不足以再容纳一个kernel A的Block,但却仍可能容纳下一个kernel B的Block.

接下来一个很重要的问题是Block如何被执行。我们可以看到,SM上的CUDA核心是有限的,它们代表了能够在物理上真正并行的线程数——软件抽象里,Block中所有的线程是并行执行的,这只是个逻辑上无懈可击的抽象,事实上我们不可能对一个任意大小的Block都给出一个同等大小的CUDA核心阵列,来真正并行的执行它们。
因而有了Warp这个概念:物理上,Block被划分成一块块分别映射到CUDA核心阵列上执行,每一块就叫做一个Warp.目前,CUDA中的Warp都是从threadIdx = 0开始,以threadIdx连续的32个线程为一组划分得到,即便最后剩下的线程不足32个,也将其作为一个Warp.CUDA kernel的配置中,我们经常把Block的size设置为32的整数倍,正是为了让它能够精确划分为整数个Warp(更深刻的原因和存储器访问性能有关,但这种情况下仍然和Warp的size脱不了干系)。
在GM204的SM结构图里我们可以看到,SM被划分成四个相同的块,每一块中有单独的Warp Scheduler,以及32个CUDA核心。Warp正是在这里被执行。
Warp的执行非常类似于SIMD. Warp中的活动线程由Warp Scheduler驱动,同步执行。我们可以看到,GM204中32个CUDA核心共享一个Warp Scheduler. 关于Warp执行中可能出现的复杂些的问题,留到下文另外说。

现在可以整理一下这个世界的图景了。SM上存活着几个Block,每个Block中的变量占据着自己的寄存器和Shared Memory,Block被划分成32个线程组成的Warp. 这样,大量的Warp生存在SM上,等待被调度到CUDA核心阵列去执行。

Warp Scheduler正如其名,是这个Warp世界里的调度者。当一个Warp执行中出现等待(存储器读写延迟等)后,Warp Scheduler就迅速切换到下一个可执行的Warp,对其发送指令直到这个Warp又一次出现等待,周而复始。这就是上一篇所说“用多线程掩盖延迟”在硬件图景下的模样。

CPU_GPU_COMPARE

图 2.  GPU用多个Warp掩盖延迟 / 与CPU计算模式的对比

本图引用自PPT “CUDA Overview” from Cliff Woolley, NVIDIA.

如图,GPU用多个Warp快速切换来掩盖延迟,而CPU用快速的寄存器来减小延迟。两者的重要区别是寄存器数目,CPU的寄存器快但少,因此Context Switch代价高;GPU寄存器多而慢,但寄存器数量保证了线程Context Switch非常快。

多少线程才能够掩盖掉常见的延迟呢?对于GPU,最常见的延迟大概要数寄存器写后读依赖,即一个局域变量被赋值后接着不久又被读取,这时候会产生大约24个时钟周期的延迟。为了掩盖掉这个延迟,我们需要至少24个Warp轮流执行,一个Warp遇到延迟后的空闲时间里执行其余23个Warp,从而保持硬件的忙碌。在Compute Capability 2.0,SM中有32个CUDA核心,平均每周期发射一条指令的情况下,我们需要24*32 = 768个线程来掩盖延迟。
保持硬件忙碌,用CUDA的术语来说,就是保持充分的Occupancy,这是CUDA程序优化的一个重要指标。

(未完待续)

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

一些后续补充。

网友邵:

SM上可以同时存在多少个Block,除了受到资源的限制之外,还受到设备上限的限制,每个SM有一个Device Limit,warps和blocks不能超过对应的上限。

CUDA, 软件抽象的幻影背后

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

今天最酷炫的事情应该就是来自老黄的这条消息:1TFLOPS,P < 15W, ARM Cortex A57 * 4 + ARM Cortex A53 * 4 +  Maxwell 256 CUDA Cores,  Tegra X1.

tegrax1
图1.  Tegra X1

本想挖掘一下写篇博,但目前报道满天飞没太大必要了。于是又想起了这个命途多舛的话题:CUDA. 关于CUDA我写了两次,第一次不满意未发,第二次成文后保存失败灰飞烟灭在热力学第二定律决定的命运里。今天借X1的东风,我们再来聊聊CUDA.

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

CUDA是个以性能为第一目标的语言,这也决定了CUDA开发者所要面对的复杂性远远要多于CUDA语言所抽象出来的编程模型本身。这大概会是软件抽象所要面对的永恒话题,我们可以去抽象出一组逻辑上漂亮完备的功能基元,却不能保证从性能的观点看它们同样也是小开销的基本操作。具体在CUDA里,最典型的例子是内存<->显存数据交换,一个简单的拷贝操作在性能上却是让人难以接受的,这背后是PCIE总线;对性能影响稍小些的例子比如Global Memory的读写需要考虑对齐,这是由于硬件层面warp和cache机制的体现;再者如过度臃肿的kernel或block过大导致寄存器耗尽,局域变量被吐到Local Memory导致的性能损失。

所有这些,都要求我们透过CUDA简洁干净的编程模型,看到软件抽象的美丽幻影背后那个不同的世界,它存在于抽象之下我们不熟悉的另一个层次,却透过性能这一个几乎是唯一的方式来影响着我们的软件。这颇类似万有引力与我们世界的关系:引力是唯一能透入额外维度的基本相互作用,如果世界有我们所不知道的维度存在,如何才能感受到那个世界对我们的影响?答案就是用引力。看过《星际穿越》的同学们想必对此有些印象。

在深入GPU的硬件架构之前,我们不妨先探讨一下这个问题:为什么GPU具有这么高的计算能力?我们试着归纳两条最主要的原因。

目前典型的计算模式有两种,CPU式的高速低延迟串行计算,和GPU式的高延迟高吞吐大规模并行计算。CPU是人们熟知的,它具有高速的内部寄存器和Cache,现代CPU又加入了多级流水线,猜测、乱序执行,超线程等技术加速其指令吞吐能力,具有快速的响应能力,但是对于大量数据的处理却相对不够用。这一点3D游戏应用就是典型的例子,当然,这就是GPU崛起的契机。
GPU天生为数据的批量处理而生,它擅长的是在大量数据上同时做同样或几乎一致(这点很重要)的计算。为什么要求一样的计算?这一点可以从很多个角度来回答。
最重要的一个回答是,多个线程同步执行一致的运算,使得我们可以用单路指令流对多个执行单元进行控制,大幅度减少了控制器的个数和系统的复杂度(设想成千上万的线程各自做不同的事情,如果再有线程间通讯/同步,将会是怎样的梦魇)。
另一方面,现实世界中应用在大规模数据上的计算,通常都涵盖在这一计算模式之中,因而考虑更复杂的模式本质上是不必要的。比如计算大气的流动,每一点的风速仅仅取决于该点邻域上的密度和压强分布;再如计算图像的卷积,每一个输出像素都仅是对应源点邻域和一个卷积核的内积。从这些例子中我们可以看到,除了各个数据单元上进行的计算是一样的,计算中数据之间的相互影响也具有某种“局域性”,一个数据单元上的计算最多需要它某个邻域上的数据。这一点意味着线程之间是弱耦合的,邻近线程之间会有一些共享数据(或者是计算结果),远距离的线程间则独立无关。
这个性质反映在CUDA里,就是Block划分的两重天地:Block内部具有Shared Memory,线程间可以共享数据、通讯和同步,Block外部则完全独立,Block间没有通讯机制,相互执行顺序不影响计算结果。这一划分使得我们既可以利用线程间通讯做一些复杂的应用和算法加速,又可以在Block的粒度上自由调度计算任务,在不同计算能力的硬件平台上自适应的调整任务安排。
现在我们把注意力放在“几乎一致”这里。最简单的并行计算方案是多路数据上同时进行完全一致的计算,即SIMD(单指令流多数据流)。这种方案是非常受限的。事实上我们可以看出,“完全一致”是不必要的。只要这些计算在大多数时候完全一致,就可以对它们做SIMD加速,而在计算分叉,各个线程不一致的特殊情况下,只需要分支内并行,分支间串行执行即可,毕竟这些只是很少出现的情况。这样,把“完全一致”这个限制稍微放松,就可以得到更广阔的应用范围和不输于SIMD的计算性能,即SIMT(单指令流多线程)的一个重要环节,这是GPU强大处理能力的第一个原因。

一个或许让每个初学者都惊讶的事实是这样一组数据:Global Memory访存延迟可以达到数百个时钟周期,即便是最快的Shared Memory和寄存器在有写后读依赖时也需要数十个时钟周期。这似乎和CUDA强大的处理能力完全相悖——如果连寄存器都这么慢,怎么会有高性能呢?难道这不会成为最大的瓶颈吗?
答案恰恰就出乎意料:不,这不是瓶颈,这个高延迟的开销被掩盖了,掩盖在大量线程之下。更清楚的说,当一组线程(同步执行,类似于SIMD的一个线程组,在CUDA里叫做warp)因为访存或其他原因出现等待时,就将其挂起,转而执行另一组线程,GPU的硬件体系允许同时有大量线程存活于GPU的SM(流多处理器)之中,控制单元在多组线程之间快速切换,从而保证资源的最大利用率——控制单元始终有指令可以发放,执行单元始终有任务可以执行,仍然可以保持最高的指令吞吐,每个单元基本都能保持充分的忙碌。
这就是GPU硬件设计中非常有特色的基本思想:用多线程掩盖延迟。这一设计区别于CPU的特点是,大量高延迟寄存器取代了少量低延迟寄存器,寄存器的数量保证了可以有大量线程同时存活,且可以在各组线程间快速切换。尽管每个线程是慢的,但庞大的线程数成就了GPU的数据吞吐能力。此为高性能的第二个原因。

这文又要写成未完待续了。接下来的日子,不填完旧坑不再开新话题。

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)这个参数矩阵使得我们可以针对各个不同的参数设置不同程度的步长约束,从而能够构造出一个稳定的算法。

(未完待补)

算法描述与性能优化的解耦——Halide语言 (1)

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

程序的结构和运行效率常常被人们看作是难以调和的。这个事实源于我们把一个数学上结构清晰良好的(比如,用递归形式刻画的)算法映射到一个现实的不完美的计算模型上,这个模型计算是有代价的,要极小化这个代价就需要尽可能的重用中间计算结果,减少依赖增加指令级并行,充分利用空间和时间的Locality和存储体系的分级结构,等等…

于是长久以来,程序中描述算法要”做什么“的逻辑,掩盖在了用来优化性能,描述”怎么做“的芜杂逻辑之下。当然,我们也有能够只通过清晰刻画”做什么“来编程的语言(如函数式语言家族,尤其是以Haskell为代表的纯函数式语言),这些语言收起了让用户自己规划计算过程的权限,把计算过程的优化交给编译器自动完成。这导致了人们对其效率的不信任,事实上这个问题也确实普遍存在,聪明的编译器只是少数,而且他们也未必能保证给出一个高度优化的程序。所以至今,对高性能有要求的各种库仍然采用C等提供了底层操作的语言实现。

性能优化会破坏软件结构这一点带来了无尽的苦恼,和风险。这使得我们每做一次优化都需要格外慎重。性能优化后的代码基本被凝固,不再具有良好的可维护、可修改性质,这样我们跨平台移植或者改变优化方案时就会面临代码需要大片重写的风险。所以人们确立了”不成熟的优化是万恶之源“这样的信条,强迫自己把性能优化留到最后一步,万不得已时采用最成熟最保守的思路去优化。

真的只能接受这个现状吗?我们知道抽象和解耦是软件的灵魂,当两件不同目的的事情纠缠在同一段代码中时,意味着需要把两件事情各自抽象出来,解耦成简单独立的逻辑。这里我们正面对一个典型的解耦问题:算法描述(做什么)和性能优化(怎么做)需要被解耦。函数式语言虽然能做到这个解耦,但是它把优化工作交给不那么靠谱的编译器了,那能否把这个优化工作交还给设计者自己呢?设想我们如果能够独立构造两个逻辑,就可以利用如纯函数式语言这样具有强大描述力的工具,寥寥数笔刻画出算法;然后再针对具体的硬件平台,实现相应的优化计算方案。不同的平台间的移植只需要替换这个优化方案部分。更好的一点是,我们可以尝试更激进的优化方案,测试各种各样的方案,这不过是个替换而已,而且对算法的功能本身没有影响。

解耦工作的难度一定程度上取决于要解耦的两个概念是否能够清晰的区分开来。算法描述和性能优化的解耦是不容易的,因为一般说来这两个概念不易区分。但在图像处理这样的领域里,计算具有典型的模式(数据在pipeline上流动,被各个节点依次处理),我们仍然可以把二者很好地解耦。

Halide就是这样一门语言。

Halide是由MIT、Adobe和Stanford等机构合作实现的图像处理语言,它的核心思想即解耦算法和优化,事实也证明它是成功的,在各种实例中它均以几分之一的代码量实现出同等或者数倍于手工C++代码的效能,更不用提代码的可维护性和开发效率。

先上例子(来自Halide的文献”Decoupling Algorithms from Schedules for Easy Optimization of Image Processing Pipelines”),三个版本的图像模糊算法,以及他们各自的性能。




void box_filter_3x3(const Image & in, Image & blury) {
	Image blurx(in.width(), in.height()); // allocate blurx array
	for (int y = 0; y &lt in.height(); y++)
	for (int x = 0; x &lt in.width(); x++)
		blurx(x, y) = (in(x - 1, y) + in(x, y) + in(x + 1, y)) / 3;
	for (int y = 0; y &lt in.height(); y++)
	for (int x = 0; x &lt in.width(); x++)
		blury(x, y) = (blurx(x, y - 1) + blurx(x, y) + blurx(x, y + 1)) / 3;
}



9.96 ms/megapixel
(quad core x86)

代码1.  C++实现图像模糊,结构良好但效率差。




void box_filter_3x3(const Image & in, Image & blury) {
	__m128ione_third = _mm_set1_epi16(21846);
#pragmaomp parallel for
	for (int yTile = 0; yTile &lt in.height(); yTile += 32) {
		__m128ia, b, c, sum, avg;
		__m128i blurx[(256 / 8)*(32 + 2)]; // allocate tile blurx array
		for (int xTile = 0; xTile &lt in.width(); xTile += 256) {
			__m128i*blurxPtr = blurx;
			for (int y = -1; y &lt 32 + 1; y++) {
				const uint16_t *inPtr = & (in[yTile + y][xTile]);
				for (int x = 0; x &lt 256; x += 8) {
					a = _mm_loadu_si128((__m128i*)(inPtr - 1));
					b = _mm_loadu_si128((__m128i*)(inPtr + 1));
					c = _mm_load_si128((__m128i*)(inPtr));
					sum = _mm_add_epi16(_mm_add_epi16(a, b), c);
					avg = _mm_mulhi_epi16(sum, one_third);
					_mm_store_si128(blurxPtr++, avg);
					inPtr += 8;
				}
			}
			blurxPtr = blurx;
			for (int y = 0; y &lt 32; y++) {
				__m128i*outPtr = (__m128i*)(& (blury[yTile + y][xTile]));
				for (int x = 0; x &lt 256; x += 8) {
					a = _mm_load_si128(blurxPtr + (2 * 256) / 8);
					b = _mm_load_si128(blurxPtr + 256 / 8);
					c = _mm_load_si128(blurxPtr++);
					sum = _mm_add_epi16(_mm_add_epi16(a, b), c);
					avg = _mm_mulhi_epi16(sum, one_third);
					_mm_store_si128(outPtr++, avg);
				}
			}
		}
	}
}



11x fasterthan a
naïve implementation
0.9 ms/megapixel
(quad core x86)

代码2.  上一段代码的优化版本,效率好但结构性破坏。




Func halide_blur(Func in) {
	Func tmp, blurred;
	Var x, y, xi, yi;
	// The algorithm
	tmp(x, y) = (in(x - 1, y) + in(x, y) + in(x + 1, y)) / 3;
	blurred(x, y) = (tmp(x, y - 1) + tmp(x, y) + tmp(x, y + 1)) / 3;
	// The schedule
	blurred.tile(x, y, xi, yi, 256, 32)
		.vectorize(xi, 8).parallel(y);
	tmp.chunk(x).vectorize(x, 8);
	return blurred;
}



0.9 ms/megapixel

代码3.  Halide代码,清晰简短,且一样高效。

我们可以看到,在Halide所实现的版本中,代码分成两部分,一部分是描述算法的algorithm部分,采用典型的函数式风格定义出要计算什么;另一部分则是指定”如何计算“的schedule部分。Halide目前没有自己的语法解析器,它的前端直接嵌入在C++里,作为一个库来使用。我们构造出一个图像处理算法后,可以把它编译到诸如x86/SSE, ARM v7/NEON, CUDA, Native Client, OpenCL各种平台上。

用schedule这样一个概念来抽象各种各样的底层优化技巧是整个方案里最关键的一环。不同平台有不同的优化技巧,怎样才能用一个统一的观点去处理它们,使得我们能在一个足够简单、与底层细节无关的世界观里处理优化问题?Halide用这样的观点来统一各种性能优化方法:它们都是控制存储或计算顺序的手段。这样,Halide通过提供一系列控制计算过程中存储和计算顺序的工具而帮助我们描绘性能优化方案。

Halide目前并没有太多的考虑编译器自动优化的问题,但这是一个漂亮的开端。如果将来在手动优化的同时仍有强大的编译器优化做后盾,将会是一番什么景象?

本站将持续跟踪这方面的进展。关于Halide语言进一步的剖析,请等待下篇。

(未完待续)

酷技术:freeD三维场景回放

昨天说了3D全景,今天再搜了下,发现了freeD这个东东。

说起来不新鲜,中文网络上这条信息也已经是一年前的了。这就是一个3D重建的典型应用,在体育场上利用多台(比如官网给出的16-28)高清相机在多个位置多个角度采集同一场景的图像,重建出3D模型。从Demo看重建质量真的不错,但不知实际运行效果如何。

视频1.  freeD三维场景回放技术

3D重建这段时间也是山雨欲来的感觉,之前放出了超炫城市3维重建Demo(视频2)的acute3D公司目前已经把爪伸到中国了,他们通过航拍视频重建出了长城的3D模型

视频2.  acute3D的巴黎三维场景重建Demo

这家公司致力于大规模的三维重建,这是个激动人心的事儿。比如目前百度地图已经有360°*90°的高清全景街景,仅用这些数据就可以重建出相当一部分城市三维面貌来,若再有航拍数据,一个真正的数字化虚拟城市也是可以期待的。如果能用微型无人机航拍采数据,可以实现廉价的大规模重建,基于此能玩出多少花样来,只是个想象力的问题。

对于想亲手玩一下的朋友,不妨试一下VisualSFM或者123D Catch,前者更学术化一些,由PBA的作者Changchang Wu开发(PBA是目前最好的开源Bundle Adjustment实现),后者是个产品。

酷技术:SamSung Project Beyond,实时3D全景

最近几个月各种实时全景拼接技术雨后春笋般冒了出来,看来一项技术到了瓜熟蒂落的时候,是挡也挡不住。今早无聊搜了下实时全景,还是把不关注技术新闻又懒于做技术推广的老夫吓了一跳。

目前市面上大多数产品跟我们类似,无非是给拼接算法一个高性能实现,或者基于FPGA,或者基于CUDA。真正让人眼睛一亮的是三星最近推出的Project Beyond,这款产品配合一个三星的虚拟现实眼镜Gear VR,可以实现真正的身临其境感——对于我们关注技术的人来说,3D,这个词儿是唯一的重点。

beyond_01

图1.  三星Project Beyond

我们知道人类之所以能感知纵深,是因为双眼上像点与场景点构成一个三角,数学上我们可以用三角测量来计算纵深,因此我们双眼感知到的信息里是包含纵深的。3D眼镜正是利用该原理,给双眼不同的图像,利用这个差异产生纵深感。

从目前的报道里,我们可以看出Project Beyond是一个赋予双目纵深感的全景装置,这一点是它超出以往技术的关键。

QQ图片20141219103205

图2.  Project Beyond的构造

如图1,Project Beyond有17个广角相机,其中1个指向天空。它采取了与普通全景相机共中心摆位不一样的摆位方式,这自然是因为要产生3D效果的要求——普通的共中心摆位无法感知到纵深信息,这可以参考我们的《全景拼接算法原理》系列文章。

 

我们也一直有把自己的实时全景技术做成微型硬件设备的想法,可惜各方面因素制约(尤其销售是我们的弱项)尚未实施。现在看起来有些可惜。

我们认为计算机视觉在接下来几年会有狂飙突进的发展——各方面条件都已成熟了,无论是理论还是硬件计算能力。而新型的人机交互手段可能是这其中最重要的一个领域,在PC迅猛发展的这二三十年里,鼠标键盘始终巍然不动,现在是时候改变了。Project Beyond这样的产品只是个开始,各种新的体验正扑面而来。

图像拼接算法原理 2

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

2. 曲面投影

Homography_Near90Degreee

图6. 近90^{\circ}时产生越来越大的畸变

通常简单的图像拼接技术,就是如上节所示的基本原理,找出一张大概处于中间位置的图像,然后利用单应性变换把其他图像变换到该中心图像的视角下,再做一些后续的曝光补偿、图像融合等处理即可。但是这一技术有相当大的局限性,最简单的例子,不能直接用它拼出360^{\circ}的全景图。

为什么呢?让我们来考虑图6中所示情形。可见,三条光线与P^{\prime}相交的三点,原本是近乎等间距均匀分布的,而当它们映射到P平面上后,间距却产生了巨大的差异。表现在图像上,P^{\prime}上的图像变换到P上后,会产生相当大的拉伸畸变。当图中两相机的投影平面越来越趋于垂直时,这个畸变越来越大,以至于P^{\prime}上普通的一点可能会被映射到P平面的无穷远点。这时,这种简单的单应性拼接方案就彻底崩溃了。

Homography_Warp

图7. 曲面投影解决大视场角下的投影畸变问题

怎样得到更宽视场角下的拼接图?解决方法很简单。以上所出现的这种畸变源自于我们将点投影到一个平面上,设想将P掰弯,或者直接弯曲成一个圆柱面,成为图7所示的样子,那原本被投影到P平面无穷远处的点就被拉回来了。我们在圆柱面上选取一个足够均匀的坐标系,把坐标对应到像素坐标,就可以得到一个全景图了。

当然,针对不同的应用,我们还可以选取不同的投影曲面,比如选取球面用于360^{\circ} * 180^{\circ}的球面全景图,甚至也可以选择一个立方体作为投影曲面。

 

3. 后续处理

至此全景拼接的几何原理就大致说完了,虽然我们还没有给出数学表达。为了先居高临下的了解整个拼接流程,我们不妨把后续处理的梗概也在此一说。

实际应用中为了创建出完美的全景图,有很多的问题需要考虑。最典型的问题有两个,一个是如何解决不同照片中曝光不一致的问题;一个是如何在拼接缝处完美平滑的融合两张图像的问题。

第一个由曝光补偿算法来解决,大体思路是估计两张图间的曝光差异,然后进行补偿。此处不多说。

第二个问题也有众多解决方案,最为著名的大概就属Multi-Band融合算法。该算法虽然八十年代就已提出,但其效果至今仍让人赞叹。在通常图像间失配程度不大的情况下,Multi-Band可以达到肉眼几乎不可分辨的融合效果。其原理也不复杂,下面略微一提。

融合两张图像,最直接的方案是在两张图像的重合区域用一个平滑渐变的权重对二者加权叠加。该方法的效果并不理想,关键原因是我们无法兼顾拼缝附近的局域细节和大尺度上两张图片的宏观特征(如光照)。当我们希望局域细节能够完好拼接时,需要用较小的平滑渐变区;而当我们希望要宏观上平滑过渡时,又想要较大的渐变区域。这二者似乎不可调和。

但事实上并非如此。Multi-Band的成功之处就是在于它同时兼顾两种需求,当融合宏观特征时,采用一个大的平滑渐变区;融合局域细节时,则采用小的平滑渐变区。那如何才能把这两种情况分开处理呢?很简单,把图像分解为不同频带的分量之加和,图像的宏观特征在它的低频分量图里,而局域特征在高频分量图里。

所以,Multi-Band算法的过程大致就是:把图像按照频率高低展开成一个金字塔,然后高低频分量各自按照不同的方式平滑加权并叠加,最后把各频带分量重新加和,得到最终的融合结果。

该算法融合效果虽好,但对于计算量要求较大,它需要创建多座金字塔并对金字塔进行各种运算,图像像素较高时,在CPU上要达到实时基本无望。当然,GPU上情况就不一样了,我们自己就实现了实时的Multi-Band融合算法,效果很好。

 

这一系列文章主要以拼接的几何原理为主。下一节开始用数学建模前两节所述的投影模型。

(未完待续)

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子空间的结构。

图像拼接算法原理 1

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明,否则将追究法律责任。
前置广告:多路视频实时全景拼接Demo可见我们的视频主页

0. 引言

76

stitch1

图1,2,3.  两张图片的拼接

图像拼接是计算机视觉中一个有趣的领域,它把来自多个不同视角相机的图像变换到同一视角下,无缝拼接成一张宽视野图像(比如360度全景图,甚至360度*180度的球面全景)。上图所示,即为两张图像的拼接,结果基本是完美的。

需要注意的是,由于相机各自的指向角度不一样,因此两图片中来自同样场景的部分并不能够通过平移图像而完全重合。比如上图中那栋带有岭南风格拱顶的房子,它的屋檐在左图中要比右图中更水平些,如果试图通过平移对齐像素来拼接两幅图,结果必然不自然(比如在屋檐处出现拐角)。两图要做到完美叠合,不是一个平移变换能做到的。当然,肯定存在这样一个变换,那它是什么呢?

事实上这里边的数学原理(射影几何)很古老,什么情况下不同位置、视角的相机可以变换到同一视角下拼接起来也是久为人知的,只不过能够利用计算机来进行大规模自动拼接,还是近些年才成熟起来的事。

那到底什么情况下图像可以拼接,如何拼接呢?不妨先摆出结论吧:在两种情况下图像可拼接,一是各相机几何中心重合;二是各相机位置任意,但场景是一个平面。一种特殊情况,场景为远景时,可以近似的等价于一个平面场景,从而也是可拼的。拼接的方法很简单,在满足上述两种情况之一时,存在一个单应性变换(Homography),能够将一个相机的图像变换到另一个相机的视角下从而可以进行拼接。

下面我们用一套直观的语言来讲解这个变换,以及更复杂些的拼接方法。

 

 1. Homography

Homography

图 4.  单应性变换

如图,计算机视觉中,我们用一个投影中心点和一张图像平面来刻画一个理想的相机。点与平面的距离为焦距,任意场景点所成像点用这样一种简单的方式来决定:将它与投影中心连线,连线与图像平面的交点即像点位置。图中给出了两个共中心的相机,它们有不同的指向(因而也有不同的图像平面PP^{\prime}),同一场景点S在两个平面上的像点由紫线与两平面的交点DD^{\prime}给出。到这里我们就直观的得到了这样一个变换:P上像点到P^{\prime}上对应像点的一个映射。这个映射可以把两张不同视角下拍摄的照片变换到同一视角下,从而能够完美拼接两张图像。这个映射正是我们要说的单应性变换。

到这里我们可以回头来解释另一个问题,为什么要求相机共中心?不妨反过来考虑,如果两相机不共中心,会出现什么样的情况?如图所示,

Homography_2

图 5.  单应性的破坏

我们让两相机的投影中心OO^{\prime}相离。这时可以看到,S_1S_2两点原本在相机O上投影到同一像点,这时在相机O^{\prime}上却投影到了不同像点。这意味着相机O^{\prime}看到了在O视角下被遮挡而看不到的内容,此时我们无法把O^{\prime}看到的图像变换到O的视角下。考虑一个极端情形可以进一步解释这个问题:如果两个相机分别拍摄同一物体的正面和背面,那它们所看到的内容无论如何也不能变换到同一视角之下,或者说,我们根本找不到这样一个新的相机,在它的视角之下可以看到OO^{\prime}二者所看到内容的总和。这就解答了我们之前的问题:非共中心放置时,各相机看到的内容太丰富以至于不能变换到同一视角下。读者可以检查,在相机共中心时,一个相机中被投影到同一像点的场景点,总会在另一相机中也被投影为同一像点。

事实上,在相机摆位任意时,我们看到的信息是如此丰富,以至于可以尝试重建出场景点的三维坐标。现代的三维重建算法可以利用不同位置、视角下拍摄的图像重建出物体的三维模型,这是后话。

当然,读者仍然可能想起,除掉共中心情形,我们还说过,相机摆位任意,但场景为一个平面时也是存在单应变换的。我们在这里不再给出这种情况的直观解释,读者可以自己思考。

 

(未完待续)