无锁编程教程:简介


无锁程序的设计非常复杂,一方面,需要读者对于并行计算涉及概念有较深的理解,另一方面,由于不同平台表现不一致,调试困难等客观原因,即使写出了程序也不容易检验。我的无锁编程始于Bruce Dawson的论文,十分幸运的是,我甚至有机会在工作中,借鉴Bruce的一些无锁数据结构编程技巧,并运用到Xbox 360等平台的编程与调试中。

从那以后,我写了一系列资料,从最基本的理论推导到最实际的代码编写,甚至最底层的硬件细节都有所涉及。然而不幸的是,这些资料甚至有些都是相互矛盾的,比如,有一些材料直接假定了顺序一致性(sequential consistency),并借此避免了任何有关于内存顺序细节的讨论,而在另一些材料中却没有做这些假设。随着C++11原子数据类型 的出现,本来就非常混乱的无锁编程变得更加复杂,让我们更难解释无锁算法的运行原理了。

在这篇文章中,我希望能够重新系统得介绍无锁编程,从对于每一个基本概念的定义开始。运用流程图,展示出这些概念如何互相关联,进而分析这些概念的一些细节。如果你准备继续读下去,我假设你已经理解最简单的多线程程序的编写,以及锁(mutex),信号量(semaphores)等同步机制。

无锁编程是什么?

大多数人对于无所编程的理解都停留在“不用锁来写多线程程序”。毋庸置疑,这样的定义是没有问题的,但是并不完整。对于无锁编程,被学术界认可的定义更加宽泛,学术界认为,“无锁(lockfree)”仅仅是一种代码的性质罢了。

简单的来说,如果程序的一部分满足接下来提到的这些条件,那么我们就可以认为他是“无锁”的,反之,如果不满足,则程序不是“无锁”的。

通过上图,我我们发现“无锁(lock-free)”这个词中的“锁(lock)”并不是指代码中的“锁(mutex)”,而是指的是能够通过某种方式“阻塞(locking up)”整个程序的可能性。这个可能性可能来源于经典的“死锁(deadlock)”,也可能由于很糟糕的内存调度顺序使得程序停滞不前(这可能变成一些恶意攻击者的突破口)。这听起来挺搞笑的,但是由于持有锁(mutex)的线程的调用顺序是任意的,因此最坏的情况是可能发生的,纵使有些时候,理论计算告诉我们“最坏情况”是一个不可能事件(概率为0),但是从概念上来说,只要“可能发生”阻塞,那么这个程序就不是无锁的。

举一个例子,对于下面的无锁程序,在初始 X=0 的情况下,什么样的调度顺序可以导致程序永远按无法退出循环?(这个问题留给读者自己思考)

没有人指望一个很大的工程是完全无锁的,通常来说,我们希望这个工程中某一个访问了公有内存的代码片段是无锁的。举一个例子,我们可能希望拥有一个“无锁队列”,这个无锁队列支持若干无锁的操作,诸如 push , pop , isEmpty 等等操作。

The Art of Multiprocessor Programming的作者Herlihy 和 Shavit倾向于将上述的操作描述成类的方法,进而对无锁(lock-free)有了一个更简单的定义:

在无限次执行中,可以完成无限次方法调用。

换句话来说,只要程序不停地调用这些“无锁”的方法,已完成的函数调用不停增加,则从算法上来看,系统是不可能发生死锁的。

从无锁编程的定义,我们可以引出一个很重要的结果:假如说我们停止了一个线程,这个操作是不会影响其他线程的执行的。正因为如此,在编写实际程序的时候,所有的线程都必须能够在有限的时间内完成,不论其他的线程干了什么。

顺带一提,如果一个操作被设计成阻塞操作,他有可能仍然是无锁的。比如无锁队列的删除会在队列为空的时候被阻塞,这并不影响他是无锁的。

无锁编程的技巧

当你尝试去实现一个无阻塞的无锁程序,这会引出一系列技巧,比如:原子操作(atomic operations),内存屏障(memory barriers), ABA问题等。从这里开始,问题就变得棘手起来。

这些技巧是如何关联起来的呢?我准备使用一个流程图来解释这个问题:

原子的“读-修改-写”操作

原子操作的定义为一系列可以访问内存,并且不可分割的操作。没有任何线程可以观测到进行了一半的原子操作。在现代的处理器中,很多的操作都是原子的,比如:基础数据类型的“对齐读取(aligned reads)”和“对齐写(aligned writes)”通常是原子的。

更近一步,我们有了 读-修改-写(Read-modify-write) (RMW)操作,这个操作使我们可以原子化更复杂的操作,特别的,当有多个写线程时,这个操作非常有用。原因在于,当多个线程试图去读-修改-写一个内存区域时,他们会自动被排成不想交的一列,并依次执行。在之前的博客中,我已经讲过读-修改-写操作了,比如在 轻量锁, a 自旋锁 和  轻量的日志系统.

RMW操作的例子有Win32编程中的  _InterlockedIncrement ,IOS编程中的 OSAtomicAdd32 ,和C++11中的 std::atomic<int>::fetch_add 。值得注意的是,C++11标准并不保证这个实现是无锁的,具体是不是通过无锁来实现依赖于你的平台。因此,你最好提前了解你的平台是否支持这些无锁操作,或者通过调用 std::atomic<>::is_lock_free 来查看。

不同的CPU系列使用不同的方式支持RMW操作。比如PowerPC和ARM等处理器使用load-link/store-conditional 指令,这种指令使你能够在底层自定义你自己的RMW操作,即使通常并没有人这么干,毕竟普通的RMW操作已经足够了。

正如流程图中所说,RMW十分重要,即使是在只有一个处理器的系统中。如果没有了原子性,一个线程可能被中途打断,并导致一些潜在的问题。

“比较并交换”循环

最常讨论的RMW操作是 比较并交换(compare-and-swap) (CAS)。在Win32里面,CAS通过一系列intrinsics函数实现的,比如 _InterlockedCompareExchange 。通常情况下,我们在循环中使用CAS操作来尝试提交一些操作。具体来说,我们首先拷贝一些公有变量到局部变量中,并进行一些操作来改动局部变量的值,最后,将这些改动借助CAS操作提交至公有可见。

这样的循环是无锁的,因为如果一个线程的if条件测试为false,这意味着他一定紧接着另一个线程执行完毕相同操作的线程之后。尽管在一些架构下,存在 弱化的CAS操作,在这些CAS操作中,上述不一定是if语句判断false的充要条件,但是这些弱化的CAS操作仍然能够有效避免ABA问题。

顺序一致性

顺序一致性(Sequential Consistency)意味着所有的线程都有共同的对于内存访问操作的顺序的判定。这个判定的内存访问操作的顺序恰好等同于在源代码中的内存访问定义的顺序。在线性一致性模型中,一些内存顺序有关的问题不会发生。

一个非常简单(但是显然不可能在实际中使用的)解决内存一致性的方式是禁用编译器优化,并且强制所有的线程在同一个处理器上面运行。这样,即使存在着线程的切换,处理器将永远无法发现自己的内存顺序出现了问题。

一些编程语言即使在优化了代码之后,在多核环境中,依然保证了顺序移植性。在C++11中,你可以定义拥有默认内存顺序(memory ordering)的公共的原子数据类型来实现内存访问的移植性。在Java中,你可以将公有变量标记为 volatile 已实现顺序移植性。下面是一个例子:

由于C++原子操作保证了顺序一致性,因此结果r1 = r2 = 0是不会发生的。从编译器层面来看,编译器添加了一些额外的指令(如内存屏障或RMW操作)来保证顺序一致性。这些额外的指令带来了更大的时间花销,但是让程序员能够以最直观的方式操作内存顺序。

内存顺序

正如流程图中所述,当你在多核计算机上运行无锁程序,且你的程序并没有保证内存顺序一致性的时候,你需要思考如何避免内存乱序的发生。

在线代的处理器结构中,这些试图矫正内存顺序的方法可以分为两类分别是矫正 编译器导致的内存乱序(compiler reordering)处理器导致的内存乱序(processor reordering) 具体而言有如下方式:

  • 一系列轻量级的同步以及屏障指令,这些我将在之后讨论
  • 一系列内存屏障的指令,我在之前已经展示过
  • 一些内存顺序的指令,提供了acquire和release语义

Acquire语义防止该语句后的语句内存乱序,而Release语义防止了该语句之前语句的错排,这些语义特别适合用于处理“生产者/消费者”关联的线程关系,一个线程生产并广播消息,而另一个线程接收消息,这些我也将在之后的文章中详述。

不同处理器有不同的内存模型

不同的处理器有不同的内存模型,在内存重拍的问题上,不同的CPU以来的规则是依赖于硬件的,比如PowerPC和ARM处理器可以修改内存写的执行顺序,而Inter和AMD的x86/64处理器则不行。因此,古老的处理器有更加弱的内存模型(relaxed memory model).

通过C++11,我们可以调用一些固有的模板来无视掉这些由于平台导致的问题。不过我认为程序员们需要知道不同平台的区别。其中比较重要的一点事,在x86/64的处理器层面,所有的读(load)操作都自带了acquire语义,而所有的写(store)操作都自带了release语义,至少对于非SSE指令集和非并行储存的内存模型而言是这样的。从结果上来看,一些在x86/64上运行正常的程序在其他处理器上崩溃就不是一件很奇怪的事情了。

如果你对于处理器如何将内存操作重排感兴趣的话,我推荐你阅读Is Parallel Programming Hard的附录C。并且时刻警醒内存重排有时候是来自于编译器而非处理器。

在这篇文章中,我们有过多的提及无锁编程的一些实际代码编写上的策略,比如什么时候我们应该用无锁编程,我们应该做到什么程度,以及验证你程序有效性的必要性。对于大多数读者,这个简介只是为了给你一些无锁编程的最基本的概念,让你在之后的阅读中不会因为一些简单的概念而困惑。与之前相似,如果你发现有任何的不准确的地方,请通过评论告诉我。

vim简单教程及配置

vim是什么?

既然今天是向大家安利vim,那么先说说安利的原因吧:

通过官网链接下载vim并且安装后,我们可以看到vim的启动界面——

唔,请关注重点……

我是在安利读者为慈善事业做贡献呢!

这个逼装的,我给99分……剩下1分哪里去了呢?

——剩下一分是因为vim确实确实非常好用

言归正传,vim是一款编辑器,而编辑器是什么呢?

编辑器就是类似于notepad(记事本)这样的只能用来打字的工具,而当这种工具集成了编译运行调试等等功能后,就成了集成开发环境(Integrated Development Environment)也就是我们常说的IDE,比如DevC++和VC6.0。

IDE与编辑器的区别就在于由于IDE集成了所有功能,所以对初学者非常友好,但是开发效率通常来说比较糟糕,反之,编辑器一心只做编辑方面的优化,所以在编写代码方面效率完爆IDE。这些在之后都会提到。

Vim作为一个编辑器,作者亲身经历表示他至少能让你的代码编写速度提升一倍,而且和vim配套使用的g++,gdb还会使你的代码调试速度变成原来的三倍甚至更多。

为什么是vim

  1. 作者作为一个vim党,当然会推荐vim咯
  2. 帮助乌干达的可怜儿童
  3. 无视个人偏好,emacs,notepad++,sublime等“编辑器”都是很不错哒(换句话来说dev-c++,vc6.0,xcode等IDE就不符合这篇文章的核心思想了),当然,vim也属于这些编辑器中的一个
  4. 以后同学们编写代码不一定是在本机,有可能需要登录到其他电脑上,这时能够在终端里面使用的编辑器据我所知也就vim和emacs了,到时候即使你是一个资深IDE党也不得不去学习使用vim。
  5. 我曾在计算机系群里的调查,结果如下(大势所趋,不服不行啊……):
    1. 明确表示使用编辑器编程的同学们中,vim党最多
    2. 对于所有的被调查者vim也有相当大的优势

vim怎么用

不同平台的vim其实都是大同小异的

我们接下来就用mac下的终端vim举例了。

安装

不同平台vim的安装方式不同

  • Windows用户:在vim官网中搜索安装方法与安装文件,需要注意的是,windows下的vim是gvim,表示图形界面版本的vim。
  • Linux用户:可以使用该发行版自带的包管理器安装vim,比如Ubuntu使用 apt install vim
  • Mac用户: brew install vim

与此同时,倘若用vim写C++程序的话,我们还需要安装g++与gdb。

  • Mac用户: brew install g++ ,gdb最好就改成mac自带的lldb好了
  • Ubuntu用户: apt-get install g++ && apt-get install gdb
  • Windows用户:
    • 方法1:首先安装Dev-C++,在安装目录中找到g++,gdb所在的目录,将目录加入windows环境变量path中,至于环境变量该怎么用,有什么意义,不是本文探讨的范畴。
    • 方法2:安装MingW,然后使用mingw-get安装g++。

运行

Windows用户

windows用户在桌面打开图标就好了。

当你看到帮助乌干达的可怜儿童那句话就说明你成功打开vim了~

注意vim存在一个当前目录的概念,windows大家最好还是通过“用Vim打开”一个目录特定文件的方式进入vim以确保当前文件夹就是你想要储存文件的文件夹。否则,倘若直接打开vim的桌面快捷方式,你保存的文件会去一个奇奇怪怪的地方。

而检验gdb/g++是否安装成功,你只需要在CMD里面输入命令 g++ ,看是否报错说没有该命令。

MAC/LINUX用户

mac、linux用户在终端输入vim命令即可,在mac/linux终端里面,你在哪个文件夹下打开的vim,vim的当前目录就是哪个。

而检验gdb/g++是否安装成功,你只需要在终端输入命令 g++ ,看是否报错说没有该命令。

vim简易教程

首先,vim官方是有教程的,在开始界面输入 :help 就可以直接进去

其次,如果想要学习vim,最好拿着电脑跟着教程亲自试一试

最后,要有心理准备,你可能会经历一周左右生不如死的适应过程。

第一步:像用记事本一样用vim

先说个笑话:

问:如何生成随机字符串
答:让萌新退出vim

大家可以自己尝试一下(windows用户不允许使用窗口那个红叉叉)就会发现,进入vim后,不但不容易退出,甚至按键盘是无法直接开始编写程序的,因为vim分为不同模式,开始时进入的模式叫做“命令模式”,是专门输入命令的,而我们需要在“插入模式”中写程序。

模式切换

命令模式:输入命令的模式。之后没有特殊指名,所以命令都在命令模式下输入

插入模式:输入字符直接插入文本的模式

  • 命令模式->插入模式: i
  • 插入模式->命令模式: <esc>

基本操作

vim中,有很多操作以冒号开头,这些操作使用频率相对较低,但不影响其重要性:

退出vim:退出 :q 或强制退出 :q!

  • 查看帮助: :help
  • 运行命令: :!<command> ,例如 :!ls 列出当前文件夹下文件(命令的执行目录为当前vim打开的文件的目录)
  • 保存文件: :w 或 :w <filename>
  • 打开文件(若不存在新建): :e filename

保存与打开文件可以为绝对路径,也可以为一个文件名,此时默认在“当前目录”打开

快捷操作

在命令模式中,命令不一定以冒号开头。

  • 剪切/删除一行: dd
  • 粘贴: p
  • 撤销: u
  • 回退: <ctrl+r>
  • 光标移动:方向键或 h / j / k / l

第二步:像IDE一样用vim

先从一个特定目录进入vim,然后编写程序,并保存为a.cpp.

我们知道,倘若在终端下切换到了源文件所在目录,执行 g++ a.cpp o a 可以将当前文件下的a.cpp编译为可执行文件,同理,调用 gdb ./a 可以使用gdb调试程序(mac是lldb)。

我们还知道,vim中,我们可以使用 :!<command> 调用外部命令,所以我们在vim中输入 :!g++ % -o %< 就等同于从终端用g++命令编译源文件。(注意:命令中的 % 表示当前vim正在编辑的文件的文件名,而 %< 则表示他去了后缀的文件名。)

通过上面所述的手段,我们可以实现在vim中的编译运行:

  • 编译: :!g++ % -o %<
  • 调试: :!gdb %<
  • 运行: :!%< 
  • 查看空间: :!size %<

vim提供了键盘映射函数map,形如 :map A : B ,那么如果我们提前用了映射命令,就不用输入编译命令了。

  • 编译: :map <F9> : !g++ % -o %< <CR>
  • 运行: :map <F6> : !%< <CR>
  • 调试: :map <F5> : !gdb %< <CR>

其中, <CR> 表示回车,读者可以试试没有这个 <CR> 会变成什么。

通过了映射,我们可以直接按 <F9> 编译,按 <F5> 调试。

最后,由于这些map命令时预处理性质的,我们可以将它写入一个叫做vimrc的文件,每次打开vim的时候自动加载

  • 在mac/linux中该文件为 ~/.vimrc
  • 在windows中,可在菜单栏-编辑-启动设置中直接修改

比如我们可以在vimrc中加入

 

第三步:定义自己的vim

之前已经知道vimrc中可以定义编译运行快捷键,实际上还可以有更多的自定义属性,比如:

  • 左边显示行号: set number
  • 查找命令高亮: set hlsearch
  • 开启自动缩进:
    • set smartindent
    • filetype plugin indent on
  • 用鼠标拖拽分栏: set mouse=a
  • 设置折叠: set fdm=marker
  • 设置缩进:
    • set tabstop=2
    • set softtabstop=2
    • set shiftwidth=4
  • 设置代码高亮: syntax on
  • 令插入模式jj为退出: imap jj <esc>
  • 设置颜色: colors evening

其中,颜色我推荐大家使用

加上这些命令之后,界面美观度一下就上升了——

最后,读者可以在文末找到我的vimrc配置文件。

第四步:神一般的编辑器

了解完之前的只是,我们也只能把vim当做记事本来用,下面介绍一些快捷键,可以极大加快编写代码的速度。在开始了解快捷键前,你需要知道vim快捷键的设计初衷是让你可以不用鼠标,甚至手不离开主键盘就可以打字。

基本操作

vim中,有很多操作以冒号开头,这些操作使用频率相对较低,但不影响其重要性:

  • 退出Vim: :q
  • 强制退出Vim: :q!
  • 查看帮助: :help
  • 运行命令: :!<command> ,例如 :!ls 列出当前文件夹下文件(命令的执行目录为当前vim打开的文件的目录)
  • 保存文件: :w 或 :w <filename>
  • 强制保存: :w!
  • 保存并退出: :wq
  • 打开文件(若不存在新建): :e filename

撤销操作

  • 撤销: u
  • 回退: <ctrl+r>

分栏

  • 竖向分栏: :vsp
  • 横向分栏: :sp
  • 不同栏跨越: :<ctrl+w> <up>/<down>/<left>/<right> 或 :<ctrl+w> i/j/k/l

多文件

打开vim时用 vim a.cpp b.cpp 或 vim *.cpp 可以同时打开多个文件,不过第一个文件全屏显示,需要通过快捷键切换。

  • 下一个文件: :n
  • 上一个文件: :N

进入插入模式

  • 当前光标前: i
  • 当前光标后: a
  • 当前行首: I
  • 当前行末: A
  • 当前行后新建一行: o
  • 当前行前新建一行: O
  • 退出插入模式:
    • 默认:  <esc>
    • 在配置文件中加入 imap jj : <esc> 之后,可用 jj 退出插入模式

进入覆盖模式

  • 进入字符覆盖模式,并在覆盖后回到命令模式: r
  • 直接进入覆盖模式: R

复制粘贴与删除

  • 复制当前行: yy
  • 剪切当前行: dd
  • 删除光标位置字符: x
  • 复制之后的10行: 10yy
  • 剪切之后的10行: 10dd
  • 全文复制: :%y
  • 复制第3行到第5行: :3,5 y
  • 剪切第3行到第5行: :3,5 d
  • 粘贴: p
  • 粘贴10份: 10p

光标跳转

  • 行首: ^
  • 行末: $
  • 第5行: 5G 或 :5
  • 光标移动:方向键或 h / j / k / l
  • 文首: H
  • 文末: L

查找替换

  • 查找: /pattern
  • 全文替换: :%s/old_pattern/new_pattern/g
  • 全文替换(每行替换第一个匹配): :%s/old_pattern/new_pattern
  • 第5至10行替换: :5,10s/old_pattern/new_pattern/g
  • 跳到下一个: n
  • 跳到上一个: N
  • 替换所有int单词(正则表达式使用): :%s/\<int\>/long long/g

代码缩进

  • 自动缩进: yy=G 或 gg=G
  • 第五行到第十行向前缩进: 5,10 <
  • 第五行到第十行向后缩进: 5,10 >
  • 接下来五行向前缩进: 5<<
  • 压行: J

其他

  • 括号匹配: %
  • 宏定义(录制操作): q<name>
  • 宏定义调用: @<name> 或 @@
  • 调用make命令: :make

再论Vim的好处

效率高

  1. 使用了vim你的手甚至可以不动鼠标,不离主键盘
  2. 编程中存在的很多类似于复制特定行,粘贴若干遍的操作在vim下都是可以用命令解决的
  3. 查找替换功能强大,连不能查找替换的很多功能都可以用宏来解决,对于批量修改帮助极大
  4. 其配套的gdb让你永远告别肉眼查错,进入二分查错的行列,复杂度直接从O(N)降到了O(logN)

用途广

  1. 基本上所有平台都有vim,而其他IDE只能在特定平台用
  2. vim的操作不依赖于具体语言,不需要每学一门语言都重新学习一个新的IDE
  3. 有丰富的插件可以使用(但是插件方面Vim比起Emacs还是差远了)
  4. Vim模式在很多地方都有支持,比如Emacs和Sublime,所以你学了Vim还可以直接迁移到这些平台上

好玩

  1. 突然发现找到了除大括号换不换行(当然要换)之外新的战场:编辑器用Emacs还是Vim(当然是Vim),引战什么的最好玩了~
  2. 我这人还是非常支持慈善事业的

结语

愿世界和平,乌干达儿童幸福快乐~

 

原文为《【和平向】一波vim安利》,计62班级公众号

配置文件

总结下来,一个简单的windows版vimrc如下(其中1-33行是自带的内容,34行之后是自定义内容),这个版本适合OI党使用,可以现场手敲vimrc。

MAC版的vimrc如下(支持了更多文件的运行快捷键)

最后,可以在这里找到我的最新的vimrc

真实感图像渲染系列:小结

这是真实感图像渲染系列的第十篇文章。

这里写一写个人在编写过程中发现的一些坑,分享出来以免之后再踩吧。

  • CMake好啊,比Makefile好,可以很容易整合第三方库,但是……编译超级慢,别人的工程3s编译完成我的需要半分钟……
  • eps设多少确实是个坑,我设的是1e-7,但是得保证牛顿迭代的精度严格高于eps,否则焦点可能存在大于eps的误差,eps的设置就没有意义了
  • MAC的错误报告挺有用的 == 比如下面这个直接秒调

真实感图像渲染系列:相机(Camera)与景深效果(Field of Depth)

这是真实感图像渲染系列的第九篇文章。

普通相机模型

不同的相机模型中,有一个人眼位置和一个屏幕,对于屏幕中每一个像素,连接人眼和该像素位置得到一条射出光线。其2D示意图如下:

首先,一直屏幕像素坐标,以及屏幕本身在场景中的位置,我们可以计算出该像素的位置,如图中红点所述。接着,连接人眼和像素的光线就可以用于计算该光线的颜色。这一部分挺容易理解的。

带景深的相机模型

在普通相机模型的屏幕和人眼位置之外,带景深的相机模型额外存在一个焦平面。如下图所示

通过普通的相机模型,我们可以求得普通相机的光线的方向。

普通相机的光线一定会和焦平面在某个位置碰撞,为了达到景深效果,一方面,我们需要保证同一个像素对应的光线一定会击打到焦平面上的同一个点,另一方面,在焦平面之外的其他位置,同一个像素对应的光线应该击打到不同的点。

解决方案如下:我们求得在普通相机模型中,光线和焦平面的焦点,将该焦点景深相机模型的焦点。接下来,我们对光线的起始位置做一个微小扰动,并且修改光线方向使得光线恒定过焦平面上的点。这样,无论初始光线如何扰动,焦平面上的点都可以最清晰的表现出来。而不在焦平面上的点,会产生随机的模糊效果。

 

 

真实感图像渲染系列:贝塞尔曲线旋转体与直线求交

这是真实感图像渲染系列的第八篇文章。

概述

贝塞尔曲线是一个由多个控制点组合而成的曲线。在本文探讨的贝塞尔曲线,是由多条简单的4阶贝塞尔曲线拼接而成。绕z轴旋转贝塞尔曲线,我们可以得到一个旋转体曲面,该曲面就是待渲染的参数曲面。本文将探讨如何使用牛顿迭代法与直线求交。求交过程中

  1. 关于迭代初值的求法
  2. 迭代中的一些优化技巧
  3. 关于浮点数精度误差eps的求法。

贝塞尔曲线的定义

贝塞尔曲线是一个由多个控制点组合而成的曲线,n阶贝塞尔曲线定义为

(1)   \begin{equation*}  P(u) = \left[ \begin{matrix}P^{(x)}(u) \\ P^{(y)}(u) \end{matrix} \right] = \left[ \begin{matrix} \sum _{i=0}^n P_i^{(x)} C(n,i) (1-u)^{n-i}u^i \\ \sum _{i=0}^n P_i^{(y)}C(n,i) (1-u)^{n-i}u^i \end{matrix} \right] \end{equation*}

式(1)定义了该曲线P(u)的关于参数u \in [0,1]的参数方程,其中P_i为n个曲线的控制点。当我看到这个公式的时候,整个人都已经懵掉了。实际上理解起来并不复杂。

函数P(u)由参数u定义了一个二维向量,其中P^{(x)}P^{(y)}相互独立。只需要理解其中一个的产生式即可。考虑函数 f_i(u) =C(n,i) u^i (1-u)^{n-i-1},发现这f_0(u),f_1(i),\cdots,f_{n-1}(u)实际上是一个和为1的数列,数列在杨辉三角形中定义。换句话来说,f_i(u)实际上可以得到一个权值的分配,使用f_i(u)对于所有的控制点P_i加权平均,就是贝塞尔曲线尝试完成的事情。当f_0(u) = C(n,0) u^0 (1-u)^{n-i} = 1权值全部分配给了P_0,因此贝塞尔曲线的起点(u=0)和P_0重合。

贝塞尔曲线旋转体

在三维空间中,我们通过旋转贝塞尔曲线得到一个曲面。曲面强制对z轴旋转对称。(注:在代码中是对x轴旋转对称,不过z轴对称相对好理解一些)。

(2)   \begin{equation*}  S(u,\theta) = \left[ \begin{matrix} Q_x + sin\theta P^{(x)}(u) \\ Q_y + cos\theta P^{(x)}(u) \\ Q_z + P^{(y)}(u) \end{matrix} \right] \end{equation*}

一个贝塞尔曲线旋转体要能够成功渲染,其核心就是需要支持与光线的求交。光线我们使用(rayO,rayD)表示,rayO表示光线起点,rayD表示光线方向,是一个单位向量。则光线上一点可表示为C(t),满足:

(3)   \begin{equation*}  C(t) = rayO + t \times rayD \end{equation*}

倘若直线C(t)S(u,\theta)相交,则有

(4)   \begin{equation*}  F(t,u,\theta) = C(t) - S(u,\theta) = 0 \end{equation*}

其中,F(t,u,\theta)t,u,\theta的函数,要求t>00\leq u \leq 1,原问题可转化为求函数零点问题。函数零点问题使用牛顿迭代解决。

(5)   \begin{equation*}  x_{i+1} = x_i - [F'(x_i)]^{-1} \cdot F(x_i) \end{equation*}

其中F'(x_i)为函数F(x_i)的Jacobian矩阵,定义为:

(6)   \begin{equation*}  \frac{\partial F}{\partial t} = \left[ \begin{matrix} \frac{\partial F_1}{\partial t} & \frac{\partial F_1}{\partial u} & \frac{\partial F_1}{\partial \theta} \\ \frac{\partial F_2}{\partial t} & \frac{\partial F_2}{\partial u} & \frac{\partial F_2}{\partial \theta} \\ \frac{\partial F_3}{\partial t} & \frac{\partial F_3}{\partial u} & \frac{\partial F_3}{\partial \theta} \\ \end{matrix} \right] \end{equation*}

要试图计算出式(\ref(eq:jacobian)),我们需要分别求出F(t,u,\theta)对于t,u,\theta偏导数的每一个分量的值(共九个)。结果为:

(7)   \begin{equation*}  \frac{\partial F}{\partial args} = \left[ \begin{matrix} rayD_x & -sin\theta \frac{dP_x}{du} & -cos\theta P_x(u)\\ rayD_y & -cos\theta \frac{dP_x}{du} & +sin\theta P_x(u)\\ rayD_z & -\frac{dP_y}{du} & 0\\ \end{matrix} \right] \end{equation*}

其中\frac{dP}{du}是对于式子(1)按照最原始的求导法则求导的结果。

贝塞尔曲面物体设计

如图

共6条4阶贝塞尔曲线,每一条由绿蓝红绿四个控制点定义。

注意事项

  1. 牛顿迭代需要保证答案中u \in [0,1],而牛顿迭代本身无法附加条件,因此需要在答案的邻域中查找合适的u,t,\theta,我用到的方法是,对于给定曲线,生成一个圆柱体包含框,随机该圆柱体里面高度h的一个圆形面片,将光线与该面片的交点的u,t,\theta作为迭代的初值。随机h约30次即可得解。
  2. 由于牛顿迭代不太精确,可能是的结果产生微小偏移,这种偏移会对判断点在平面哪一侧产生影响,我们可以通过调整eps的取值,不同部分赋予不同的eps,使得这些误差不至于相互影响。
  3. 贝塞尔曲线求交本身费时间,但是判断是否一定没有交点却相对简单。可以求出贝塞尔曲线旋转体的“圆柱体包围盒”,(更精确的话,可以是空心圆柱体包围盒)。然后判断直线是否与包围盒有交点。判断直线与圆柱交点P = rayO + k \times rayD的位置,可以先讲所有东西映射到z=0屏幕,这是变成2维直线和圆求交,算出k,回带到三维情况验证。
  4. 实际上,像本文采用的多条低阶贝塞尔曲线求交并不明智,倘若将本文那个碗的6条贝塞尔曲线改为1条简单一些的曲线,然后做一个玻璃杯什么的,渲染效率将提高6倍!

 

 

真实感图像渲染系列:使用OpenMP并行优化渲染

这是真实感图像渲染系列的第七篇文章。

如果将渲染修改为并行,那么渲染效率将提升约8倍(我的电脑是八核)。非常非常有诱惑力。而且基于OpenMP框架的并行化甚至还非常简单!

平台编译

使用OpenMP最麻烦的地方是框架安装问题,倘若本机的g++版本支持 -fopenmp ,那么恭喜你,你已经完成了90%的工作了。倘若本机g++不支持该编译选项,那么需要在网上搜索对应的安装方法。安装方法不同系统都有较大的差别,故不在这里细讲。

OpenMP程序并行化指令都是用 #pragma 实现的,也就是说倘若编译器不支持,也不会产生什么问题。

OpenMP简单使用

这是一个最简单的for语句,对应于PPM算法中的光子映射。

该for语句满足如下条件:

  • 循环起始位置和终结位置确定
  • 不同的循环没有依赖性,

满足这样条件的for语句,理论上,将循环不同的i并行处理是没有问题的。OpenMP使用了最简单的预编译指令帮助我们实现for语句的并行,修改为

OpenMP临界区设置

#pragma omp critical 编译器将自动将 #pragma 指令后面紧接着的语句块进行并行。

一般来说,并行语句块内使用的变量,都是读写局部变量以及读全局变量,这时,不会存在什么问题,但是总有特殊情况,我们多个并行的循环需要修改同一个全局变量。这是可能导致错误的。如PPM算法中,计算“视点”,需要将结果保存在全局数组 view_pts 里面。这时,多个线程并行修改会导致一些问题(诸如SegmentationFault)。这时,使用 #pragma omp critical 可以强制后面紧接着的那个语句块通过“加锁”保证同时只有一个线程在运行该代码块。

 

真实感图像渲染系列:渲染可视化

这是真实感图像渲染系列的第六篇文章。

动机

一幅图片渲染时间少则一分钟,多则半天,等到渲染完成再看效果未免也太佛系了。况且诸如PPM算法亮度调节什么的一般在前两轮的时候,就已经大概知道了,并不需要渲染完。因此我们希望写一个实时显示渲染效果的类。

实现

显示部分在技术上没有任何难度,就是OpenCV中的imshow而已,如果必要,每次imshow之前,使用imwrite备份图片到文件。

我们需要一边显示一遍渲染,这个就需要多线程了。由于opencv的限制,imshow和waitKey只能在主线程调用,因此我们把图像渲染放在子线程里面去。这里线程的创建我使用了pthreads库。

具体细节就看代码啦。init函数传入一个Color类数组的指针,而PaintBoard能够实时监测该指针指向内容的变化,并且显示出来。(需要用户在窗口按回车键刷新)

 

 

真实感图像渲染系列:渐进式光子映射算法(Progressive Photon Mapping)

这是真实感图像渲染系列的第五篇文章。提前声明,这篇文章讲述的算法之和PPM在思想上相似,细节上由于笔者也没有仔细研究论文,所以可能会有比较大的出入。

算法概述

光线追踪算法 中,我们提到了光线追踪实际上无法成功捕获到一部分光路。随之而来的是透明物体下方出现错误的阴影。

渐进式光子映射算法(Progressive Photon Mapping),尝试分别从摄像机和光源两边一起投放射线。两边的射线在某漫反射面上相聚。RayTracing和PPM的关系挺像深度优先搜索和双向深度优先搜索的关系。

摄像机出发有若干视线,这些视线落在漫反射面上形成视点。由于视线本身存在归属的像素,所以每一个视点也存在归属像素。假设该视点位置得到了光强,那么该光强将累加至对应像素。

与此同时,光源处随机发出若干光子,光子经过折射,反射,以及随机漫反射后,也会在若干位置碰撞漫反射界面。我们希望光子会对碰撞位置周围半径R的所有视点有一定的光强贡献。称之为“送光”。

从摄像机出发视点的计算,和RayTracing视线的计算相似,唯一的区别是RayTracing在视线碰撞漫反射面后,通过采样光源得到此处的光强,而PPM在得到的所有视点后,使用一颗KD树进行储存便于后续的“送光”。

从光源处发的光子,其起始位置、方向皆随机选取。与视线相比,唯一的不同是光线并不会终止于一个漫反射面,在漫反射面上,光子仍然会随机选取方向反射,并且,在对于该漫反射面碰撞点周围R距离内的所有视点,都进行“送光”。

对于整体稳定的考虑

PPM算法中,光子是随机发射的,由于光线追踪本身速度并不快,所以有限时间内,实际可以发射光子数是有限的。另一方面,如果我们场景是一个“点光源”和一个“地平面”,很容易发现,点光源随机发射的可以达到地面可见区域的光子更少。我们现在需要考虑的,就是如何合理利用这些足够少的光子,使得渲染出来的图片尽量均匀。

下图是R比较小导致送光不均匀的例子,观察地面上“斑驳”的光影:

光子本身需要给半径R周围的所有视点送光,一方面,提高R有助于使得送光范围更大,少数离散的光子也能渲染出均匀的图片。另一方面,将送光光强定义为半径R从中心到两边线性递减的函数,可以消除边界处光强骤减的问题。

对于局部高亮的考虑

过大的R导致局部细节无法突出,下图是R较大的情形,可以发现橙色透明玻璃球的聚焦并不算成功。

因此最优的方法应该是将渲染分成若干不同的round,R随着渲染的进行而逐渐变小。这样,前期图像趋向于亮度均匀,后期,细节被强调出来。

对于整体光强的考虑

每个光子究竟需要携带多少光强这是一个非常难以衡量的问题。由之前的讨论,我们意识到了光子的发射需要区分round,不同round对应的R不同,换句话来说,前面的round和后面的round同等重要。但是,倘若前面round和后面round光子拥有相同的光强,这会下面的问题:

由于像素本身颜色是对应所有视点接收到的光强总和,在RayTracing算法中,我们限制光子的“折射”“反射”皆存在系数,最终的光强一定是小于等于(1,1,1)的。在PhotonMapping由于是随机化算法,我们是无法保证这一点的,沿用RayTracing的光强算法,我们得到的像素光强将是一个无限接近于(0,0,0)的值,这并不是我们想要的。同时,如果将所有像素光强乘以一个值,也会导致问题,比如一块正常的蓝色(0.1,0.1,1),如果不幸被乘以了系数10,最终得到的是(1,1,10) = (1,1,1),蓝色就变成白色了。这个问题可能导致下图所示的结果,不是我们想要的。

解决思路很巧妙,核心是调整每一轮的光强,我是用的方式是将每一轮光强设置为调和级数。即最终图片光强为1 \times round_1 + \frac{1}{2} \times round_2 + \frac{1}{3} \times round_3 + \dots + \frac{1}{n} \times round_n. 实际上这是一个巧妙地想法,因为后期的round期望对整体亮度不产生改变,而对局部亮度改变,恰好调和级数之和不收敛,但是单纯取出奇数或偶数项的话,他是收敛的假装他就突然收敛了。意味着倘若一个位置在之后的每一次都有很多光子碰撞,则该位置倾向于变量,而倘若之后光子碰撞频率较低的话,该位置亮度倾向于收敛。

实际上,究竟多少光强是一件很灵活的事,例如我发现高亮部分不明显,那么我可以定义每一轮初始光子光强随轮数增强,增强的系数等同于半径衰减系数等等。即R' = R \times \alphaE' = E\times \frac{1}{\alpha}.

KD树的建立

我工程中使用的KD树是基于维度轮换分割的,写起来非常方便。这里单独讲一讲KD树编写中的一些技巧:

  • KD树建树期间,按照指定维度分割出较小的一半和较大的一半。使用std::sort函数排序固然可以,但是存在更加高效的算法:std::nth_elements(begin,begin+k,end)可以提取序列第k大,将小于第k大数的放在前面,其余放在后面。将k设为n/2,则KD树的分割部分直接通过调用nth_elements完成。
  • KD树通过对子树维护点集包围盒,实现快速判断一个询问点是否在包围盒中。实际上,很多人即使得到了包围盒以及一个点的坐标,仍然无法得到一个高效的算法判断点距离包围盒的距离。这里介绍一个比较快捷的方式:
    • 令包围盒用maxX,minX,maxY,minY,maxZ,minZ定义,询问点(X,Y,Z)
    • dx = max(0,|X-maxX|,|X-minX|),dy = max(0,|Y-maxY|,|Y-minY|),dz = max(0,|Z-maxZ|,|Z-minZ|)
    • 答案dist满足:dist = \sqrt{dx^2+dy^2+dz^2}

渲染效果图

源代码

progressive_photon_mapping.h

progressive_photon_mapping.cpp

 

真实感图像渲染系列:光线追踪算法(Ray Tracing)及基于Hash的抗锯齿方法

这是真实感图像渲染系列的第四篇文章。

算法概述

光线追踪算法可以算是图像渲染中的Hello World了,不过这个Hello World可没有那么好写。通过前两篇文章,我们已经知道了渲染算法包含的核心要素,如颜色,物体等表示。

算法本身网上已经有很多详细的教程,多说无益,本质上RayTracing就是最直接的通过视线反推光路,并忽视一切不容易追踪的光路的算法。

光线追踪本身可以处理的光路需要满足两个限制:

  1. 只经过了一次漫反射面
  2. 漫反射面是作为摄像机到光源光路中最后一个界面存在

我们可以通过一个玻璃球在地上是否会聚焦光线形成光斑来判断一个算法是不是普通的RayTracing。这里,漫反射面(地面)并没有作为光路的最后一个面(中间还插入了玻璃球的投射面)。光线追踪算法会在玻璃球下产生一个黑色的阴影。

算法步骤

下面,我将以最简练的方式概述RayTracing的算法。

对于屏幕上的每一个像素,我们可以通过摄像机(camera)找到一条视线光路,RayTracing算法目的就是得到这条光路的颜色。当然,倘若光路直接射到了光源内部(即与光源产生碰撞),那么颜色就是光源颜色。否则,倘若光路射到了无穷远处,那么我们将提前设置好的背景光颜色作为实现颜色。

多数情况,视线既不会与光源碰撞,也不会到达无穷远处。视线沿着这个光路行进,在物体表面可能由于“折射”,“反射”等而分叉,我们递归得处理折射,反射后视线的颜色。这里,递归总深度,光线强度皆可作为停止递归的信号。

最终,部分视线将到达一个漫反射平面,这时,光线亮度设置为能够直接照亮该位置的光源亮度。对于复杂光源,诸如面光源,我们需要在光源处随机取点,计算出光源中多少光可以不受遮挡到达反射面。

基于hash的边缘超采样

通过基本的光线追踪算法,我们可以得到如下的图片

可以发现,球体的边沿存在较严重的锯齿效应。通过基于hash的边缘超采样后,效果如下

锯齿处的颜色变味了两边颜色的混合,从而在肉眼看来消除了锯齿。

在RayTracing中,实现的函数调用可以看做是一棵树结构,每一个节点都会碰撞一个物体。在每一次碰撞后,我们可以维护光路hash值hash' = hash \times 19 + object_{hash}。通过如下方法迭代出来的hash值,只要中途光路的折射顺序有任何细微差别,都会导致hash值的变化。我们对于摄像机每一个像素的hash值进行比较,倘若一个像素四周hash值都和他相同,则其不在边缘上。否则,我们对该像素(x,y)进行超采样。其颜色为(x,y),(x+1/3,y+1/3),(x-1/3,y+1/3),…,(x-1/3,y-1/3)这九个像素位置对应颜色的平均值。

在实际运行中,超采样的运行时间相对会比较长。毕竟稍微复杂一点的场景,会有很多边缘点,而每个边缘点的重采样时间是原来的9倍。

源代码

ray_tracing.h

ray_tracing.cpp

 

真实感图像渲染系列:使用json初始化场景

这是真实感图像渲染系列的第三篇文章。

场景初始化的必要性

不同于普通的编程作业,图像渲染的场景结构非常复杂,比如同样是物体,一个球体需要读入中心坐标以及半径,而平面则只需要一个定位坐标以及法线方向就能够确定。如何初始化一个场景,为每一个参数赋予正确的初始值是非常重要的。

在参考过往代码发现场景的初始化无非两种方式:一种是将场景作为一个类写在源代码中;另一种是将场景写在文本文件中,对于所有需要初始化的类写一个接口来读取文件。在我的代码中,使用了第二种初始化方式。就编译效率上来说,第二种方式可以防止每次修改模型都重新编译。这在后期可以优化出很多时间。

场景初始化方式

一方面,我使用json格式储存场景,json格式的读取在C++中确实不算简单,需要安装jsoncpp库。jsoncpp库提供了一个数据类型 Json::Value ,这个数据类型支持整数下表访问(用于访问列表),以及字符串下表访问(用于访问字典),以及各种转换函数如 val.asDouble() 等。最终场景的json文件附在本文末尾。

jsoncpp的简单使用可以参考 通过jsoncpp库的CharReaderBuilder解析字符串

我使用伪·visitor设计模式初始化场景,对于每一个需要初始化的类申明一个成员函数 void accept(const Json::Value& val); 接受一个 Json::Value 用于初始化自身。

例如,对于vector来说,accept函数如下。

而一个包含vector的sphere类,可以拥有如下accept函数

可见,至少在安装了库之后,json的读取是足够简单的。

为场景添加功能

场景本身会被后续不同的渲染算法所使用,为了增加代码复用性,我们可以把渲染算法公共的一些函数放在场景中,例如求一个光线最近碰撞到的物体的函数 findCollidedObject 等。最终场景类如下

后续的光线追踪等代码,都会直接继承scene来使用scene的这些接口。

源代码

场景json文件