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

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

动机

一幅图片渲染时间少则一分钟,多则半天,等到渲染完成再看效果未免也太佛系了。况且诸如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文件

 

真实感图像渲染系列:颜色,向量与物体

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

颜色

颜色,向量与物体是各种渲染算法的基础。

颜色定义为三个取值在0-1之间的浮点数,分别表示RGB分量。由于分量的值恰好为0-1之间的浮点数,颜色不仅可以理解为每种颜色的强度,还可以理解为每种颜色光线的穿透率。换句话来说,这样定义的颜色同时支持乘法和加法。

特别的,对于透明物体的折射,光线穿过物体,存在一定的吸收,该吸收率可以看做是距离的指数函数。即Color_{absorb} = e^{-dist \times Material_{absorb}},如在红光吸收率0.1的物体中穿过2单位距离,则射出的红光应该为原来射入红光强度的e^{-2*0.1}.

向量

向量的定义,也是三个0-1之间的浮点数,表示xyz坐标分量。向量同时具有表征方向和表征位置的功能,坐标本身在某种意义上也是坐标原点O加上某个表示方向的向量。

一个光线,我们就可以用两个向量rayOrayD分别表示他的起始位置以及方向,rayD是一个单位向量,某点P在光线上当且仅当P = rayO + k*rayD,且k>=0

 

物体

由于C++的多态机制,物体本身是一个抽象类,只需要支持:

  • 判断和某个光线的交点
  • 判断某个位置对应的贴图颜色
  • 提取出诸如漫反射率等表面材质的性质
  • 询问物体hash值(用于抗锯齿等)

碰撞

碰撞是渲染里面非常核心也很复杂的一个模块,我们使用一个专门的类记录碰撞的信息,实际上,对于碰撞,可以通过碰撞位置C,碰撞法线方向N和如何光线方向I完全表述。有了这三个量,我们可以简单的计算出这个碰撞的反射、折射方向(其中折射方向需要物体折射率)。这里讲一个小的注意点,就是反射,折射光线的起始点不应该都设置为C,反射需要保证光线的起始位置和入社光线位置同侧,而折射恰恰相反,由于精度误差,可能C实际在物体表面的某一侧而我们用C作为新光线的起点是不妥的。解决方案很简单,对于反射,我们使用C_{surface} = C+N*eps,折射则使用C_{backface}=C-N*eps作为光线的其实位置。

 

 

真实感图像渲染系列:工程概述

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

这是图形学课程的一次大作业,基础要求是实现光线追踪算法(Ray Tracing),包括反射,折射以及漫反射,以及参数曲面求交(Bazier Curves),在此基础上,我增加了渐进式光子映射算法(Progressive Photon Mapping)以及景深效果(Depth of Focus). 不过由于本人毫无美术细胞,所以场景设计一团糟,只是粗略体现出了每一个算法的现象而已。结果如图

包含两个含贴图的球体,一个看起来像水果盘的贝塞尔曲线,一个透视的玻璃球,以及一个距离摄像头太近的,对焦失败的球。

工程源代码可以在 https://github.com/mhy12345/RayTracing-ProgressivePhotonMaping 找到。

其中,该工程依赖了glog(实际上master版本不依赖,但是需要他进行一些必要的调试信息输出),eigen(用于矩阵运算),opencv(用户实时显示渲染结果以及bmp格式保存)以及jsoncpp(用于解析json文件)。后面三个库都需要提前安装在电脑中。同时,工程需要你的编译器支持 -fopenmp 选项,例如,在我的笔记本上,编译器 clang-omp 支持该选项,编译命令即为

根目录下的 ./serial 和 ./main 分别是渲染程序的串行版本和并行版本。

python音频处理库librosa教程(2) hop_length的选取

教程中很坑的一点,是让我们自定义hop_length,然后……what is hop_length?

现在,我们需要重新思考一下特征提取的密度,对于一个22050HZ的采样数据,显然,最后提取出来的特征序列不能比22050HZ还要密集。当最后特征等于22050Hz时,提取出来的特征当然就是自身,当特征频率更小,比如2205Hz,那么我们就可以将连续的10个数据进行做平均值、方差等等操作,使得特征长度缩小至1/10。

这种分窗口压缩数据的方式确实很容易想到,通过将连续的若干个帧数据进行一些操作,合并成一个向量,即可表示这一段时间每一帧的特征。而这个帧数变少的序列就是frame。

因为音频采样数据本身就是在一个波形上的一些随机点,这些随机点进行平均值这类传统操作操作没有任何意义!怎么办呢,其实结果很简单,既然我们希望对于一个窗口提取音频信号,那么我们直接对于这个窗口进行一次傅里叶变换就好了。因此这也就是librosa所做的,对于大部分样本提取函数,你都需要传入一个hop_length,作为傅里叶变换的长度,显然,而这个值一般选择2的整数次幂(为了可以通过快速傅里叶算法优化性能)。

对于不想想太多的我们,完全可以直接使用每个函数的默认hop_length进行编写。而教程中人为设置hop_length确实是非常有误导性,让人以为这是一个必须自己设置的量。

python音频处理库librosa教程(1) 一个简单的音频特征提取程序

这是librosa教程的第一篇。主要通过解释官网的样例程序来解释如何使用librosa提取音频特征。librosa是一个音频特征提取库,在论文librosa: Audio and Music Signal Analysis in Python中提出的,如果在学术中使用的话,可以引用该论文。


本文的目标是在python中实现对于音频特征提取,提取出一个时序的特征向量。

librosa是一个专用于处理音频的库,其官方教程可以在 Librosa官方教程 中找到,这篇教程也是基于官方教程写出来的。

笔者学习一个库一般都是直接看一篇样例程序,一般来说,弄懂样例程序究竟讲了什么,也就弄懂了这个库了。作为一个目标为提取音频特征的代码库,官方索性就直接用一个简单的“音频特征提取”代码作为样例程序。这也给我们带来了很大的便利。

音频读入

首先,我们需要知道声音是如何在电脑中储存的,声音信号作为一个连续的函数,在计算机中,只能够用采样的形式储存。

具体来说,采样就是将原来的连续函数中,指定间隔求出一系列离散的时刻对应的信号大小,如librosa默认每秒钟采样22050个采样点,则原来的音频信号就变成了每秒钟22050个值得离散函数。显然,当采样点足够密集,那么我们得到的序列就能够足够逼近原函数。在波形图中,声波平均振幅可用来表示响度的大小,而每一个周期的时间(过零率)可以用来表示音高。

librosa.load 函数的具体实现细节如下

返回值的第一部分,就是将连续的音频信号按照给定的采样率(sample rate) 进行采样的结果,这个结果是一个一维的 numpy 数组,可用于后续的分析。

返回值的第二部分,就是采样率,也就是说完全是函数参数中sr的拷贝。在少数情况下,我们设置函数参数为 sr=None ,采样率就完全依照MP3默认的采样率返回,这个时候返回值 sr 就是我们评判数组 y 代表的音频究竟有多快的唯一方式了。

由简单数学推导可以得到音频时间T可以由ysr决定

    \[T = \frac{len(y)}{sr}\]

波形分割

接下来,y作为一个合成波形,可以分成两个分量,即谐波(harmonic)与冲击波(percussive)。由于笔者也不太清楚他们具体该怎么翻译,所以按照自己的理解自由发挥咯。

从粒度上来看,谐波相对为细粒度特征,而冲击波为粗粒度特征。诸如敲鼓类似的声音拥有可辨别的强弱变化,归为冲击波。而吹笛子这种人耳本身无法辨别的特征,是谐波。接下来,我们将分别对于冲击波以及谐波进行特征的提取。感性理解,冲击波表示声音节奏与情感,谐波表示声音的音色。

节拍提取

beat_frames指的就是每一拍(beat)所对应的帧(frame)位置。如果不明白beat,frame这些术语都是指的什么,那么可以参考官方文档.

这一句话真的是一个非常有误解性质的语句,对于hop_length的详细讲述可以在 python音频处理库librosa教程(2) hop_length的选取 里面看到。倘若无法理解的话,这里给出一个捷径:对于librosa里面的所有函数,都会有一个默认的 hop_length=512 ,也就是说,我们完全可以删掉这句话,对于之后的hop_length,我们都直接使用函数的默认值就好了。

节拍提取在音频特征提取中的意义在于,节拍是所有曲子共有的,不同节拍(beat)中,音频特征或多或少有些不同。在我们模型中,如果我们每隔1ms提取出了一个音频特征(如音高),我们不如将同一拍的所有音高取平均值作为这一拍的音频特征。这样表征效果会好很多。

音频特征提取

接下来,我们可以计算梅尔频率倒谱系数(MFCC),简单来说,MFCC可以用以表达曲目的音高和响度的关联。经典的MFCC的输出向量维数是13,也就是说在非必要的情况下,这个 n_mfcc 参数就不要改了(这是笔者投paper的时候用血换来的教训啊)

MFCC本身只反映了音频的静态特性,所以我们需要对他进行差分,以得到动态特性。即音频是“如何”变化的。

既然我们已经分别得到了mfcc以及mfcc_delta,那么我们可以将它们合称为一个特征,合并过程只是单纯的二维矩阵拼接,这里用了np.vstack函数。接着,librosa.util.sync函数用于以窗口形式合并多个连续的变量,比如这里,就以beat_frames为分界点,对每一个分界点中的序列用其平均值表示。到了这一步,我们得到的是,长度为音乐拍数,宽度为特征数的矩阵。

下面是对于谐波进行分析

由于这里出来的是一个类似于色度图的东西,所以笔者非常怀疑其在后续训练模型的意义是否足够。之后,用同样的方法,将色度图按照beat_frames进行简化。

附录

节拍提取算法详情:Beat Tracking by Dynamic Programming

基于Librosa的拓展:实时的音频特征分析库rcaudio