标签归档：计算机视觉

玩玩三维重建

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

我们在实时三维重建方面的工作今年已经密集展开。或许不久后某一天，你会在本站看到带有SLAM（即时定位与地图构建）功能的四轴飞行器，或者让你在书桌上打一场现代战争的增强现实应用。在敲锣打鼓欢天喜地亮出我们自己的三维重建实现前，先拿别人的东西给大家打打牙祭。

中科大刘利刚教授的3D建模软件与处理软件简介介绍了N多实用的3D相关软件。而基于照片的快速建模软件并不多，之前玩过123D Catch，很赞。围着你要建模的物体拍摄一圈，用123D Catch加载拍摄的图像，经过其强大的处理能力，生成具有纹理的3D模型。下图是我重建的我的蒙奇奇。你要做的只是拍照、上传、等待而已，相当简单。

但是123D Catch也存在一些局限，完全的黑盒子，对重建过程没有任何操控力。

要想了解从照片如何一步步重建出三维模型，并能操控某些过程，可用的免费开源软件也不少，较常用的是VisualSFM和Meshlab：

第一步：VisualSFM

VisualSFM软件允许我们上传一系列图像，它从这些图像中找到每一个图像的特定特征，利用这些特征信息重建出3D模型的稀疏点云，而后还可进行稠密点云重建。

输入：围着要重建对象拍摄的一系列照片；

输出：一个 .out文件，存储着每个相机的位置及重建出的稀疏点云；

一个.ply文件，存储着由稀疏点云重建出的稠密点云。

第二步：Meshlab

可用Meshlab对3D网格/点云做各种操作。输入VisualSFM的生成文件，Meshlab通过一系列操作可创建出包含纹理的、干净的、高分辨率的网格，并自动计算UV映射及创建纹理图像。

输入：VisualSFM的生成文件，.out文件和list.txt文件（存储照片序列）；以及.ply文件；

输出：一个.obj文件，3D模型的网格；

一个.png文件，任意大小的纹理图；

完整的流程见下图：

第一步：运行VisualSFM

1. 输入一系列图片

拍照注意事项：切忌不要站在原地，仅转动身体去拍：相机共中心能拼接全景，但是给不出三维重建的深度信息。要以待重建的对象为中心，围着它每转10-20度拍一张，这样转一圈，有不同高度信息更好。VisualSFM没有照片数量限制，照片越多，重建出的细节越丰富，但重建过程花费时间越长。

2. 特征检测及匹配

因照片可能存在旋转、缩放或亮度变化，此过程利用SIFT算法提取、描述特征，用 RANSAC算法过滤掉误匹配。此过程亦可利用GPU加速。工作状态实时显示在侧边的log窗口。

QQ图片20150314232955 QQ图片20150314233141

3. 利用SFM进行稀疏3D重建

利用 SFM 方法，通过迭代求解出相机参数和三维点坐标。即重建出3D模型的稀疏点云。若有“bad”相机（位置错误或朝向错误），结合工具栏上的“3+”按钮和手型按钮即可删除之，使结果更准确。

4. 利用 CMVS/PMVS 进行稠密3D重建

CMVS/PMVS需自己下载，编译，也可直接下载exe文件。而后把pmvs2.exe/cmvs.exe/genOption.exe文件放到VisualSFM.exe的同目录下。

通过 CMVS 对照片进行聚类，以减少稠密重建数据量，而后利用PMVS从3D模型的稀疏点云开始，在局部光度一致性和全局可见性约束下，经过匹配、扩散、过滤生成带真实颜色的稠密点云。(下图为用Meshlab查看效果图)

至此，VisualSFM的工作告一段落，结果都已存盘。若因图片匹配失败或图片较少导致某区域重建失败或重建出的某区域细节不足，可以返回添加一些这个区域的照片，重新来过（本人较懒，未作补充，谅解）。但因特征检测和匹配的结果已存盘（每张图像对应的.sift 和 .mat文件），所以已经匹配好的图像不必再次匹配，会更快完成。

第二步：运行Meshlab

1. 打开bundle.rd.out 文件

a. 按钮1，打开由 VisualSFM生成的存储在xx.nvm.cmvs文件夹下的 bundle.rd.out 文件。随后会询问是否选择照片列表文件，选择同文件夹下的 “list.txt”即可。这一步会把相机及对应的照片导入进来，对后续的纹理处理至关重要。

b. 按钮2，打开显示层目录，检测相机载入是否正确， Render –> Show Camera，因可视化相机的尺寸比网格尺寸大得多，所以需调整相机的缩放因子，scale factor可以从0.001开始调小，直到相机位置清晰可见。

2. 稠密点云代替稀疏点云

a. 按钮3，隐藏可视的稀疏点云；

b. File –> Import Mesh加载稠密点云（xx/00/models/option-0000.ply）；VisualSFM生成多个.ply文件时，需合并成一个mesh。在载入的任何一个.ply上右键选“Flatter Visible Layers”。

3. 清除杂点

按钮4选中杂点区，按钮5删除之。大致清了桌前的一些杂点。

4. 网格化

Filter –> Point Set–> Surface Reconstruction: Poisson.

利用Poisson Surface Reconstruction算法由稠密点云生成多边形网格表面。

参数可调， Octree Depth：控制着网格的细节，此值越大细节越丰富但占内存越大运行起来慢，一般设10，可慢慢调大。

Poisson表面重建算法会生成一个“不漏水”气泡，把所有场景对象包裹在其中。即模型是封闭的。可以移除多余的面Filters –> Selection –> Select faces with edges longer than，而后删除。

保存（整个project和mesh）。

5. 修复流形边缘

后续的纹理处理要求网格化的模型必须是流形(MANIFOLD)的，因此需删除非流形边（简单讲就是任何由多面共享的边）。Filters –> Selection –> Select Non-Manifold edges，而后删除之。

6. 参数化（Parameterization）

Filter –> Texture –> Parameterization from registered rasters。

根据相机投影关系创建UV映射。

QQ图片20150314235004

保存 (整个project和mesh)。

7. 投影纹理

Filter –> Texture –> Project active rasters color to current mesh, filling the texture。

可设置任意分辨率（512的2的二次方倍:512 / 1024 / 2048 / 4096 / 8192…）的纹理图。

QQ图片20150314235126

6和7可以合为一步： Filter –> Texturing –> Parameterization + texturing from registered rasters.

8. 完成、导出

当你调整满意了之后，File –> Save mesh as… a .obj文件。即可便有了一个包含你选定分辨率纹理的obj文件。

QQ图片20150314225720 QQ图片20150314225733

收官啦。而后关乎应用，就是拼想象的时候了！

更多细节参见：We Did Stuff 。

TransProse：将经典名著转换成音乐

发表评论

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

这两天《平凡的世界》电视剧上演，N多人甚至习大大都点赞，又一股重温经典小说风。《平凡的世界》好像是我高中课堂窝在桌洞里看完的。名著曾是苦逼学生时代的调味剂。老哥被奴役的时候还编程分析过不同世界名著的语言风格。今天看到一好玩的事儿，TransProse项目把不同的世界名著转换成对应情感的音乐，听了一下，真心不错。

TransProse 读取小说文本，通过文本分析确定八种不同的情绪（快乐，悲伤，愤怒，厌恶，期待，惊喜，信任和恐惧）和两种不同的状态（正或负）在整个小说中出现的密度。音乐同步小说，按时间顺序分为beginning, early middle, late middle, and end 四部分，每一部分都有对应的音乐表示：利用情感密度数据根据不同的规则和参数来确定音乐的速度，调，音符，八度等。详见其paper。

算法描述与性能优化的解耦——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三维场景回放

1条回复

昨天说了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实现），后者是个产品。

图像拼接算法原理 2

1条回复

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

2. 曲面投影

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

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

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

图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_i$ 中 $d_0$ , $d_1$ … $d_{i-1}$ 方向的分量）。原因在于 $r_i$ 自然的和 $D_{i-1}$ A-正交，因此构造 $d_i$ 只需要从 $r_{i}$ 中扣除掉 $d_{i-1}$ 分量即可。

因此，共轭梯度算法有两处关键，一处为通过利用A-正交关系来把搜索问题转换到一个特殊坐标系，使得各搜索方向上的极小值问题解耦；一是利用 $r_i$ 和 $D_{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_i$ 由 $r_i$ 和 $d_{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_i$ 在 $D_i$ 的补空间里。这点是明显的，因为每一个优化步都会优化掉 $e_j$ 在相应 $d_j$ 方向上的分量，因而剩余的 $e_i$ 就只位于 $D_i$ 的补空间里。于是我们有

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

这样 $e$ 和 $d$ 之间的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

2条回复

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

前置广告：多路视频实时全景拼接Demo可见我们的视频主页。

0. 引言

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

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

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

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

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

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

1. Homography

图 4. 单应性变换

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

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

图 5. 单应性的破坏

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

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

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

（未完待续）

Mind Spectrum

Formerly Singularity Vision

标签归档：计算机视觉

玩玩三维重建

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

第一步：运行VisualSFM

1. 输入一系列图片

2. 特征检测及匹配

3. 利用SFM进行稀疏3D重建

4. 利用 CMVS/PMVS 进行稠密3D重建

第二步：运行Meshlab

1. 打开bundle.rd.out 文件

2. 稠密点云代替稀疏点云

3. 清除杂点

4. 网格化

5. 修复流形边缘

6. 参数化（Parameterization）

7. 投影纹理

8. 完成、导出

TransProse：将经典名著转换成音乐

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

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

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

酷技术：freeD三维场景回放

图像拼接算法原理 2

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

2. 曲面投影

3. 后续处理

Notes on Conjugate Gradient Method

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

线性方程组和极小化问题

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

共轭梯度法的理解

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

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

图像拼接算法原理 1

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

前置广告：多路视频实时全景拼接Demo可见我们的视频主页。

0. 引言

1. Homography

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

第一步：运行VisualSFM

1. 输入一系列图片

2. 特征检测及匹配

3. 利用SFM进行稀疏3D重建

4. 利用 CMVS/PMVS 进行稠密3D重建

第二步：运行Meshlab

1. 打开bundle.rd.out 文件

2. 稠密点云代替稀疏点云

3. 清除杂点

4. 网格化

5. 修复流形边缘

6. 参数化（Parameterization）

7. 投影纹理

8. 完成、导出

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

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

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

2. 曲面投影

3. 后续处理

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

线性方程组和极小化问题

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

共轭梯度法的理解

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

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

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

前置广告：多路视频实时全景拼接Demo可见我们的视频主页。

0. 引言

1. Homography

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