Epipolar a personal journal

Ceres Solver 的编译

Ceres Solver 的编译比起 SuiteSparse 就更加需要耐心了。主要原因是它的依赖比较复杂。其中有两个要点:

  1. Ceres Solver 依赖 SuiteSparse ,也就需要一个 BLAS/LAPACK 。
  2. Ceres Solver 需要 glog ,否则性能会打折扣。

先说第二点,这个比较简单,从官网下载 glog 的代码,然后编译。我使用的是 Visual Studio 2015,于是 glog 的编译就不开心了。在向 Windows 移植的过程中,glog 自己实现了一些用来做 I/O 的函数,比如 snprintf ,VS2015 提供了兼容的函数,这就产生了冲突。解决的方法也很简单,把 glog 自带的禁掉。或者从 github 的项目中下载最新的 commit,其中提供了 CMakeList 可供直接编译。

然后来解决折磨人的小妖精 BLAS 。Windows + BLAS + CMake 永远是个难题,是个难题,是个难题(重要的问题说三遍)。原因有这么几个:

  • BLAS 实现多又多,在 Windows 下好用的没几个
  • CMake 的 findBLAS 和 findLAPACK 各种不好用
  • 文档匮乏

关于 BLAS ,我使用了 Intel MKL,它的最新版本已经可以免费下载,并且在 “Wintel” 平台下性能基本上是最优的。因为 Ceres Solver 自己采用了 OpenMP 等 API 进行了并行处理,如果底层的 BLAS 也独自进行并行,可能会产生错误的行为,所以在 Ceres Solver 中使用 Intel MKL 必须要注意使用它的串行版本。

为了使用这一版本的 MKL,要到 Intel® Math Kernel Library Link Line Advisor 上获得应当链接的正确的库,参考配置如下图:

与此同时,如果使用的是 MKL 的动态链接版本,日后使用 Ceres 的程序总是要带着各种 DLL 的小尾巴,显得很不方便。结合我自己的情况,我选择使用静态链接的版本,上图中对应了 Select dynamic or static linking: Static 。

接下来,来到 Ceres Solver 的编译,在 CMake (cmake-gui)中打开代码树并初次进行配置。接下来就会见到大量的提示,按部就班地配置好 glog 和 Eigen 的路径,我们开始怼 CMake 的 BLAS 配置。

首先,为了要让它出现 BLAS 的配置,我们需要告诉它“我要用 SuiteSparse!”,接下来 CMake 就会去找 SuiteSparse,然后发现它依赖 BLAS,然后发现找不到,然后来傻乎乎问你 BLAS 在哪儿,从而诱导出 BLAS 相关的配置。

……但是故事远没有这么简单,一般情况下,你会发现 CMake 提示了各种 NOT FOUND……

冷静,这个时候就需要一些小小的魔法,怎么做呢?在 CMake 的界面里点击那个 Add Entry,添加一个叫做 BLA_STATIC 的布尔值,设置为 True 。

这样便告诉了 CMake:我要使用静态版本的 BLAS 。加上这个之后,再不要忘记打开 SuiteSparse 的开关(之前没找到 BLAS ,于是 CMake 就自作主张给你关了……),重新配置。

……然后你会发现,旧的 NOT FOUND 还在,还出现了三个新的 NOT FOUND ……

不过冷静对比一下,你会发现新出现的 NOT FOUND 恰好对应了我们在 Link Line Advisor 中得到的三个文件名。把这三个文件的完整路径对照填进去,再次打开 SuiteSparse ,配置。

经过这么一番折腾,终于见到了曙光。如果一切正确的话,CMake 便会提示找到 BLAS 和 LAPACK,并为你打开 SuiteSparse 的相关配置(又是一大堆的路径等待你人肉)。全部添加好,再一次配置,如果全部通过,就可以生成用于编译的工程了。

关于如何在 CMake 里优雅的使用 MKL ,还可以参考我之前的一篇文章 CMake & Intel MKL 之黑魔法

Visual Studio 下编译 SuiteSparse

工程中要使用来自 Google 的 Ceres Solver 。它可选依赖于 SuiteSparse ,加入之后可以大幅提高大型稀疏系统求解速度。然而在 Windows + Visual Studio 中编译 SuiteSparse 并不是一个容易的事情。之前我进行了一次编译,这里记录一下,供日后参考。

SuiteSparse 是使用 C 语言编写,本身采用了 Makefile 进行编译。为了在 Windows 下进行编译,一种办法是安装一个带有 GNU Make 的环境,例如 Cygwin 。通过分析 Makefile 的内容,可以发现它基本上是将每个源代码文件编译后链接成一个库。是非常标准的做法,因此我们可以用 Visual Studio 来完成这个过程。但是 Visual Studio 自带的 C/C++ 工程并不支持 SuiteSparse 代码中用到的一种特殊技巧——基于宏的 C 代码模板。

在 SuiteSparse 中,经常出现下面的特殊用法:

#if defined(LONG)
#define Integer long long
#else
#define Integer int
#endif

此后,通过在编译时命令行加入宏定义来控制编译得到使用 long longint 的两份二进制。这种技巧需要对同一份源代码编译多次。 Visual Studio 的 C/C++ 工程不支持多次编译,也就不能直接做这个事情。

解决这个问题的方法之一是为需要多次编译的代码创建副本,从而每个副本编译一次。为了避免大量的文件复制,可以在工程中创建下面这样的源代码文件:

#define LONG
#include "suitesparse_source.c"

当然建立工程时一定要检查每个 SuiteSparse 自带的 Makefile ,确认哪些源代码应用了宏模板,哪些只进行一次编译。

针对 SuiteSparse 4.5.3 版本,我用上面的方法建立了 Visual Studio 2015 的工程,请见 SuiteSparse-MSVC 。利用它,可以不借助 GNU Make 完成 SuiteSparse 的编译。也就减少了部分配置环境的负担。

当然,这不是长久之计,更应该对 SuiteSparse 的项目管理工具进行升级,改用诸如 CMake 等的工具自动根据需要生成上面的包装代码文件。

(封面图来自 Prof. Tim Davis

从平面引申出的射影变换(二)

最近各种事情比较忙,一直没有更新……接下来陆续恢复更新。

很久以前记过一篇从平面得到两视图间射影变换的文章。在那篇文章中,平面在第一个相机坐标系下的法线是 $n$,到第一个相机中心的距离是 $d$,第一个相机到第二个相机的坐标变换为 $R$、$T$。那么我们知道这个平面上的点从第一个相机的射影坐标到第二个相机的射影坐标可以由射影变换 $H=R+\frac{1}{d}Tn^T$ 得到。

考虑一个特殊平面 $z=0$ ,这个平面到第一个相机中心的距离是 0 。此时不存在从第一个相机到第二个相机的射影变换。

在这个平面上建立 X-Y 坐标系与世界坐标系的 X-Y 重合,这个平面上的(二维)坐标到第二个相机投影平面上的坐标存在一个射影变换:

将 $R$ 按列记作 $(r_1, r_2, r_3)$,就可以得到:

也就是说,从这一平面上的二维坐标到相机 $R,T$ 上的射影坐标之间存在一个射影变换 $H = (r_1, r_2, T)$ 。

从内参矩阵到 OpenGL 投影矩阵的转换

做立体视觉时经常要将结果绘制出来,这时可能会面临这样一个问题:相机的内参如何转换为 OpenGL 中的投影矩阵?

http://www.songho.ca/opengl/gl_projectionmatrix.html 有一个非常好的 OpenGL 投影矩阵的介绍,不过如果跟随这个介绍来尝试转换内参则显得有些繁琐。这里介绍一种相对简单直接的方法。

我们知道:对于三维空间中的点 $(x,y,z)$ ,内参矩阵 $K$ 可以将它投影到图像像素坐标系下的一个齐次坐标。我们把它记在这里:

而投影后的像素坐标对应于 $(x_p, y_p) = (x’/z, y’/z)$

OpenGL 中,屏幕空间采用的是一个标准化设备坐标系(NDC),详情参考上面的介绍。这个坐标系并不是像素坐标系,因此我们需要进行适当的转换。

对于横纵坐标,这个转换非常的直接,假设视口的(像素)大小为 $w\times h$,那么我们就有 $x_{ndc} = 2x_p/w-1$,$y_{ndc}=2y_p/h-1$ 。将这一关系在齐次坐标下表达,就可以得到:

在 OpenGL 里,所有的坐标都是三维齐次坐标,绘制时前三个分量会除以第四个分量。将上面的关系对应到三维齐次坐标下的变换矩阵,就对应有:

注意到,这里空出了一行,这一行对应了 $z_{ndc}$ 。而它的处理与前两个分量有些不同……


在 OpenGL 的 NDC 坐标系下,深度范围是 [-1, 1] 。其中 -1 对应了一个最近的深度 $n$,1 对应了一个最远的深度 $f$,这两个深度边界都需要人工指定。

那么假如这两个边界给定了,我们只要线性地将 $z$ 映射到 [-1, 1] 不就好了吗?错!我们来看看为什么不能这么做:

我们在齐次坐标系下进行线性的映射,那么就有 $zz_{ndc} = Az+B$ 。将两边同时除掉 $z$ 就会得到 $z_{ndc} = A + B\frac1z$ 。这表明 $z_{ndc}$ 是关于 $\frac1z$ 线性的,而不是关于 $z$ 线性。将上面的两个边界条件代入求解 $A$ 和 $B$ ,可以得到:

或者

所以,完整的 OpenGL 投影关系是


由于深度上存在这样一个倒数的线性关系,如果最远深度过大,在最远深度附近,$\frac1z$ 的变化可能会非常小。这时如果深度精度不足,这种小变化可能由于舍入误差而产生错误,导致错误的深度结果。这种现象被叫做 z-fighting 。在上面的链接中有更详细的分析。

注意会发生这种现象是因为我们采用了线性的齐次坐标变换,在这种变换下,深度的映射必然包含了倒数关系。在传统的静态 OpenGL 管线中,我们只有设置投影矩阵来进行三维点向视口的投影,因此无法避免这种可能产生 z-fighting 的行为。但现在多采用可编程管线了,因此可以采用下面的方式进行投影,使得深度关系是线性的:

uniform mat3 K; // intrinsic
uniform vec2 S; // viewport size
uniform float n; // near
uniform float f; // far
varying vec3 vertex; // vertex coordinate
void main() {
    vec3 p = K*vertex;
    vec3 result;
    result.xy = 2*(p.xy/p.z/S) - 1;
    result.z = 2*(p.z-n)/(f-n) - 1; // linearly map depth
    gl_Position = vec4(result, 1);
}

注意的是上面只是演示一下思路,代码还有优化的空间。


细心的话,可能会注意到一个隐藏的小细节:上面我们得到的线性变换矩阵的行列式是大于零的,这意味着如果我们原始的坐标处于一个右手系,那么变换后的坐标也应该是一个右手系坐标。然而 OpenGL 的 NDC 采用的是一个左手系,这样我们变换后的坐标套用在这个坐标系下就会发生一个镜像变换。所以我们应该在变换矩阵中考虑这个镜像的因素,保证变换后的 NDC 坐标不会发生镜像。

我们来看一看这个镜像出现在什么地方:对于像素坐标,我们一般采用左上角为原点,向右下为 XY 正方向的朝向,如果我们建立相机局部坐标系时也采用同样的 XY 朝向建立右手系,可以得到 X 向右,Y 向下,Z 向前。而在 NDC 中,Y 是朝上的,这一差异带来了镜像问题。因此我们需要将变换矩阵中 Y 对应的行翻转一下(取负)。

不过如果是做离屏绘制,并且在绘制后将像素缓冲区/纹理缓冲区的内容读取到外部,这个镜像又无需考虑了。这是由于 OpenGL 缓冲区中,像素是逐行自下向上表示的,而外部的像素缓冲通常是自上向下表示的。因此如果将完整的缓冲读取出来,并直接套用外部的表示,图像内容将发生上下翻转。这一巧合刚好弥补了前面提到的 NDC 坐标系手向差异带来的镜像问题。因此读取出的图像又不存在翻转了。

补充:注意,因为存在提到的镜像问题,在绘制时,面的朝向会发生反转。如果原先采用 GL_CCW (默认如此)为正向,现在应当改为 GL_CW

OpenGL 的 Texture 和 RenderBuffer

在 OpenGL 中使用 Framebuffer 进行绘制时,需要为它指定绘制到的目标。可以用的目标有两类,分别是 Texture 和 RenderBuffer 。为什么会有这两种目标来做同样的事情呢? 只有经过细心地比较才会发现:OpenGL 中的 Texture 只支持色彩和深度缓冲区格式,并不支持如 Stencil 等辅助的缓冲区格式。RenderBuffer 则支持了所有的缓冲区格式。 此外,在建立两种对象的时候,需要为它们指定数据类型。对于 Texture ,需要指明数据的内部逻辑类型和一个具体的存储格式(比如24位深度 GL_DEPTH_COMPONENT24);对于 RenderBuffer 只需要指定内部的逻辑类型就可以了(比如同样储存深度,只需要用 GL_DEPTH_COMPONENT)。这说明 RenderBuffer 是一种更偏向图形硬件内部表达的对象。 在使用上,两者除了对应的 API 有细微的区别,其他方面基本可以相互替代,但 Texture 可以作为绘制目标允许我们进行 Render to Texture 。

后记:后来发现,根据内部类型/用途的不同,RenderBuffer 也可以指定具体的存储格式。