Low poly in CUDA

什么是Low-poly?

最近在和队友一起做Parallel大作业,课题是用Low-poly风格化的并行实现。说到Low-poly,就是把一张图片变成各种单一色彩组成的样子。在做课题的时候,我按照了Wenli Zhang[1]的paper来实现基本的步骤,分为三个,选点三角化光栅化。每一个步骤分别先实现了线性的版本,然后放到CUDA中去加速。没有用到OpenGL这种自带GPU加速的框架。最后的效果图如下。因为时间的关系并没有优化边缘的部分,所以会有黑边。

选点

选点,顾名思义就是在图片中选择一系列的点,用作下一步三角化的顶点。这其中会遇到两个问题:

根据Zhang的Paper,维护图片的原有信息可以先对图片实现边缘检测,再用不同的权重分别选取边缘和图片上的任意点。这样可以保证图片内部的原有形状可以被比较好的保存下来。因为对边缘检测的要求不是很高,所以就只用了最简单的Sobel算子来做边缘检测。

Sobel算子的实现很简单,首先把图片灰度化,然后就是一个简单的图像卷积。Sobel的公式如下。对于每个像素,他的值就是周围8个像素与他自己应用如下公式的结果。

对于Sobel算子这种图像卷积,我们非常容易就可以将他实现到CUDA中来优化。因为每一个像素的工作都是独立的,不宜赖于任何别的操作。数据之间也不会有联系,竞争。所以就只要把这部用像素为单位并行出去就好。gridDim就是ROUND_UP(w/blockDim.x, h/blockDim.y).

因为并行的级别是像素,Sobel通过CUDA可以实现非常高的加速度,如图所示。

三角化

三角化是整个项目中最难的部分。三角化的目标是在空间中把点连起来。但是我们也对于这些三角形也有一些基本的要求

为了实现这点,我们需要用到Delaunay Triangulation(DT, 德劳内三角化),或者 Constrianed Delaunay Triangulation(CDT)。 因为DT有两个非常重要的性质

而CDT是保证大多数三角形满足DT条件,极少数三角形可能不保证DT的性质。对于Low-ploy渲染这种非严格数学建模来说,应该是可以接受的。

常用三角化算法

传统的DT算法有插值法和分治法。插值法非常容易理解,有许多库比如delaunay.js主要分成

  1. 生成一个能覆盖所有三角形的超三角形,然后将这个三角形加入candidate_list
  2. 一次插入每一个点,对于每一个点
    1. 找出每一个违反性质的三角形,将他们移除candidate_list,将他们的边加入open_edge_list
    2. 去除open_edge_list中所有重复的边
    3. open_edge_list中的边与插入点构成三角形,加入candidate_list
  3. 删除任何包含超三角形顶点的三角形。

而分治法的思路是先将点做成有序,再将整个平面不停分割,直到每个区域只有3个点包括以下。连接每个区域的点,再将不同的区域merge在一起。因为这个涉及到一些计算几何的知识,所以我们并没有实现这个算法。

然而这些算法放到CUDA中并不能很好的发挥,主要是因为:

  1. 插值法的依赖太多,每个点都依赖上个插入的结果。可以并行的区域只有验证三角形那一块。如果并行插入到各个区域,又会出现多个worker访问一个三角形这样的数据竞争问题。
  2. condition的语句太多,不合适发挥CUDA最好的性能,因为CUDA是data parallel,SIMD的模型。Branch太多会导致严重的Divergent,从而不能发挥SIMD最好的性能。
  3. 有竞争就要加锁或者原子操作,也会影响性能。

尝试了用OpenMP将插值法并行化,得到了效果不是很好,在很多点的情况下,才可以得到5.4x的加速。在点数量不够多的情况下,启动线程的花费甚至会使运行时间更长。

CUDA 三角化

在CUDA上做三角化实现了Meng Qi[2]一篇算法中提到的方法。使用他的方法可以很快的形成CDT。 主要的思路就是先找到生成Voronoi图,因为voronoi图在数学上被证明是DT的对偶问题。

然后再找到图中的链接点,最后用prefix sum来计算每一个线程应该更新的buffer offset,做到并行地在buffer中插入三角形。其实paper还有好几步像是edge flip, 生成convex hull这样步骤。在这种project中我们仅仅粗略地实现了最初的几个步骤,得到的效果还可以就没有完成所有的步骤。最后的并行加速是非常明显的。

Voronoi

在paper中,生成Voronoi图使用了Jump Flooding算法,概念就是每个点一次向外传播自身的信息,当空像素收到别的点的信息时,他会判断这个点是否是离他最近的。如果是的话,会将自己也染成那个点的信息,并且帮助那个点向外传播消息。

这一步可以在CUDA中如果以点为单位进行并行,就会出现数据竞争,因为多个点可能尝试去更新同一个像素。然而如果以像素级别来并行,让每个像素查看自己周围是否有点,并且决定染哪个颜色,则可以避免这种问题。所以gridDim就是ROUND_UP(w/blockDim.x, h/blockDim.y).

数三角

当染完色生成voronoi图时,我们就可以通过遍历每个像素以及和他相邻的像素来判断这个像素是否是三个以上点的连接处。这一步可以通过简单的数颜色来实现。如果有三个颜色,则说明这个像素连接了三个点,则能构成一个三角形,如果是四个颜色,则说明有四个点,可以构成两个三角形。在这一步我们可以统计出每一行有多少个三角形,并行的单位是行,所以gridDim就是ROUND_UP(y/blockDim.x).

连接三角

当我们知道每一行有多少个三角形之后,就可以使用prefix sum来算出每一行应该在哪个buffer offset进行插入。之后我们就可以以行来做并行,对于每一行,找到所有的连接点,然后通过读取相邻像素的信息,就可以找到哪几个点需要被相连,然后直接构成三角形。

对于CUDA三角化的这几个步骤,还可以使用jump-list等来进一步减少遍历的时间。

光栅化

最后的光栅化是要让每个在三角形内部的像素都被渲染成三角形重心处的颜色,这一步在CUDA中的实现也非常简单,因为数据之间的依赖几乎没有。主要有两个思路

  1. 以三角形为单位来并行,找到每个在三角内的像素,并且将他们渲染。
  2. 以像素为单位来并行,覆盖该像素的三角形,并且将他们渲染。

因为像素的数量比三角形多的多,而GPU又适合data parallel,所以以像素为单位来并行会有更好的效果。以三角形为单位来并行还要考虑到边缘像素的数据竞争。 以像素为单位的并行最简单的方法是每个像素遍历所有的三角形。但对于三角形数量比较多的情况下,可以将图片划分成好几个区域,先并行的计算出落在这个区域内的三角形有哪些,然后对于区域内的每个像素,只需要遍历这几个可能的三角形就好。使用这种方法能让循环的次数变少很多。

总结

在这个项目中学到了很多CUDA编程的技巧,并且还直观的了解到什么算法合适并行,什么算法不合适。怎么样的并行粒度是比较好的。总是是一次非常有意思的体验!

Reference:

[1] Zhang W, Xiao S, Shi X. Low-poly style image and video processing[C]//Systems, Signals and Image Processing (IWSSIP), 2015 International Conference on. IEEE, 2015: 97-100.

[2] Qi M, Cao T T, Tan T S. Computing 2D constrained Delaunay triangulation using the GPU[J]. IEEE transactions on visualization and computer graphics, 2013, 19(5): 736-748.

Ivan Yang 03 May 2018
blog comments powered by Disqus