阅读：图像风格化三角剖分

Stylized Image Triangulation - COMPUTER GRAPHICS FOURM

Kai Lawonn 1 and Tobias Günther
link:https://doi.org/10.1111/cgf.13526

总览

论文提供了一种图像三角剖分的算法。相比之前的方法，该方法在考虑匹配图像特征的同时考虑到减小风格化后的图像和原图像之间的标准差，即将标准差看为目标能量函数，最小化这个函数，这样就变成了一个优化问题¹。

为了缩小标准差，论文使用一种算三角形顶点和标准差之间的梯度的方法去不断调整顶点，直到给定的标准差要求被达到（即，标准差小于某个阈值）。由于算梯度是个很费时的操作，因此论文还提供了一种用GPU实现的在像素之间并行的（rasterization based）算法来快速算梯度。

同时，为了保证艺术家能控制算法的运行结果，论文还提供了直接/间接的工具让你能插入和删除三角剖分顶点，反转三角边，或者对某个多边形区域进行重新三角剖分。

最后，论文选了几种三角剖分的渲染方式（比如常数颜色填充，也就是一个三角形填一种颜色，或者第顶点插值，即三角形三个顶点每个顶点有一种颜色，三角形内的颜色在三个顶点间插值插出来，等等）来展示最终结果。论文还试着去风格化视频，来展望这个三角剖分的潜力。

前人的工作

现有的大量开源与商业化软件里的风格化三角剖分算法大多基于：

1. 找到特征点，特征点可能是物体边缘之类的表示图像元素特征的东西。
2. 对特征点进行三角剖分。

举个例子，一些算法会先对图像过一个模糊算子（也就是所谓的低通滤波），避免特征点密集地产生在图像的一处噪音太强或者频率过高的地方，然后跑某些算法生成特征点，最后跑三角剖分。

这样的方法往往让图像的柔和且色调统一的区域被巨大的三角覆盖，而边缘以及对比度强烈的区域会生成大量的小三角形。很多时候用户还能通过一些参数调整特征的密度，从而间接控制三角剖分中三角形的数量。

然后还有一些其他风格化三角剖分的算法和相关的数值算法，论文也有提到，还提了很大一托，因此我当然没有去看，所以这里就不说了。然后这个论文表示，我很棒棒，因为我是第一个将标准差作为“能量函数”，把三角剖分看成一个“最小化目标能量函数”的优化问题去做的。

定义能量

首先定义图像函数I:U\subset R^2\rightarrow C，在这里$C$不是复数域，而是[0, 1]^3 \subset R^3或者单灰度的计算机色彩空间，而$U$是图像全集（一般来说就是个长方形框框）。

定义三角剖分的三角形几何是\mathcal{T}，满足\bigcup_{T_i\in \mathcal{T}}=U，其中T_i是三角形（三角形覆盖的点集）。同时，对于任意T_i, T_j\in \mathcal{T}，满足T_i\cap T_j是且仅是一个顶点或者一条三角形边。

对于一个三角形T，定义三角形颜色函数f_T: T\rightarrow C，f_T(x, y)可以是一个线性函数（在三角形渲染使用顶点颜色插值时可以把f_T表示成一个线性函数），在最简单的情况下，也可以是一个常数函数f_T(x, y)=c，此时每个三角形有且只有一种颜色。

写了那么多，终于可以定义我们要最小化的能量了，也就是标准差，然而论文上文说时标准差，到这里却变成了平方差，不过没关系，写出来：E(T) = \int_T (I(\textbf{x})-f_T(\textbf{x}))^2 dA

由于我们无法直接通过一些参数控制E(T)，所以我们再强行对每个三角形定义一个能量附加p(T)，，这个p(T)可以包含一些自制的参数来帮助我们（以及艺术家们）控制最终效果。同时定义一个这个附加的权重\lambda（越大，那么风格化算法越受制于我们自己的参数）。先别管p(T)具体是什么，首先提出算法的最小化目标：

E=\sum_{T\in \mathcal{T}} (E(T)+\lambda p(T)) \tag{1}

梯度下降

好，能量用数学形式表示了，接下来考虑怎么最小化这个能量。发现这玩意既然能用数学形式表示出来，而且好像形式很简单，那么自然是可以算梯度的。最近因为机器学习大红大紫的 梯度下降 方法就可以直接移植过来。

首先观察到每个顶点只影响其相邻三角形的能量，则可以把$(1)$改写成顶点形式的。

\begin{array}.
E=\sum_{\textbf v\in\mathcal{V}} E(\textbf v) \\
E(\textbf v) = \frac{1}{3}\sum_{T\in \mathcal{A}(\textbf v)} (E(T)+\lambda p(T))
\end{array}\tag{2}

然后画风一转，我们突然回过头来使用人类智慧考虑$p$中应该塞什么因子的问题。我们的算法可以通过对每个顶点计算顶点位置相对于能量的梯度，然后把顶点向能量降低的方向移动，从而期望得到更低的能量。但这样会产生一个问题：顶点可能移动到非法的地方。考虑一个顶点能移动到的位置：所有和这个顶点相连的顶点组成了一个多边形。为了方便思考，我们假设²它们组成了一个凸包。那么引入一个被称为拉普拉斯平滑（Laplacian smoothing）的修正因子，塞到p(\textbf v)里：

p(\textbf v) = [2\vert \mathcal N_1(\textbf v) \vert]^{-1}\sum_{\textbf w \in \mathcal N_1(\textbf v)}\vert\textbf w - \textbf v\vert^2 \tag{3}

这个和吊打GTY那道题目比较像，最小化这个的点应该是多边形的某个中心，应该也是一个二维的单谷函数，即，只有一个极值点。应该是有严格数学证明的不过我没有查。其中\mathcal N_1(\textbf v)被称为\textbf v的1-环（1-ring），应该在拓扑学上有定义的，但没去深究。

然后就是喜闻乐见的梯度下降了。定义一个步长h，那么对于每个三角剖分中顶点\textbf v的更新公式是：

\textbf v = \textbf v+h\frac{\mathrm{d}E(\textbf v)}{\mathrm d\textbf v} \tag{4}

h当然可以像ML里一样随时间衰减，由于这里已经开始套梯度下降的板子了，所以和梯度下降有关的优化应该也可以用进来，但这个我不懂，如果咱们做这个题目，这个点上应该能搞点结合与创新。

同时注意到每个顶点是独立的，因此可以放到GPU上去并行。

初始三角剖分

论文提到了两种生成初始三角剖分的方法，一种是直接均匀地方格撒点：

另一种是根据某个概率密度函数撒点，然后跑Delaunay三角剖分，这个概率密度函数可以通过用户界面画出来，或者通过前人研究出的特征提取算子搞出来。

优化

优化即论文中的Refinement³ 。在上文的讨论中，梯度下降过程中三角形顶点总数是不变的，但为了达到更好的效果，我们接下来将使用自适应进化打破这个规则。

对一个三角形T，当(\frac{E(T)}{\vert U \vert})^{\frac{1}{2}} ⁴大于某个用户设定的阈值时，在三角形的质心插一个新的顶点。delaunay三角剖分的一些构造算法本身就是通过不断插入顶点运算的，自然能支持在插入这个顶点后进行一些调整使得整个三角剖分满足（或者接近）delaunay三角剖分的要求。

容易注意到梯度下降过程会破坏三角剖分的delaunay性质，但这没办法解决——毕竟在符合delaunay三角剖分与达成能量最小这两个条件之间一定要有一个取舍。论文的做法是——二者交互进行，一波梯度下降，一波进化加上调整，再一波梯度下降...

注意这里提到了一个新操作”调整“，调整指删除一些面积太小，乃至为负数⁵的三角形（因为梯度下降可能会产生这样的三角形），具体操作是 坍缩 调三角形最短的一条边，将这条边链接的两个相邻顶点合并为一个顶点。

其它

后文基本就是讲用户界面和交互以及实现细节了（比如梯度具体怎么算的，让三角形染色函数$f$弄成线性函数的扩展方法...），如果我们选这一个论文去实现可能得看。后面还稍稍介绍了一下在视频风格化的应用：由于这个算法本身就是基于梯度下降不断调整的，天然适合对连续影像进行风格化，从而生成视觉上连续的画面。然后稍微分析了一下这个算法的收敛率，这个没细看。最后是性能分析，实时是不可能的，怎么都不可能的，实际上论文的代码实现在512x512的lena上配合一块1080每次迭代都要20毫秒。但是就最后效果来看还是很不错的。

估计是Artist编辑三角形编辑的不耐烦了，才会有这么粗糙的作品...不过就但论结果这个论文做的是很不错了。

总结

效果的确很棒，论文有它的演示视频大家可以去看。算分项目准备实现一个简单的这个东西，CUDA不太妙（只有N卡支持，且在不同环境上跑很麻烦），于是GPU交互准备用OpenGL4，现代的图形接口都支持Compute Shader之类的东西了，应该是可以的（但还得去学一波）。

补充

上面提到过三角形颜色填充可以有很多方式，最简单的就是常数填充，要最小化方差肯定填平均值。但是如果能够插值填充（即f_T(x, y) = ax+by+c），为了确定a, b, c，就得摄入一点点数学了。

我们想要最小化E(T)，使用像素去近似，令三角形T覆盖了q个像素，像素集合为S，那么

E(T) \mathop{=}^\sim \sum_{\textbf x \in S} [f_T(\textbf x)-I(\textbf x)]^2 \tag{5}

假设这个函数拥有一切好的性质（连续，可导，等等），极值点在导数为的地方取到，然后根据直觉极值点一定是最小值，于是算偏导数：

\frac{\partial E(T)}{\partial a} = \sum_{\textbf x\in S}\frac{\partial[(xa+yb+c-I(\textbf x))^2]}{\partial a} = 2\sum_{\textbf x\in S} [x^2a+xyb+xc-xI(\textbf x)] = 0 \tag{6}

类似的可以对b, c分别求偏导数，这里略去，整理之后最终获得一个3x3的线性方程组，可以表示为：

\left[\sum _{\textbf x \in S} \left( \begin{matrix}x^2 & xy & x \\ xy & y^2 & y \\ x & y & 1 \end{matrix} \right)\right] \times \left(\begin{matrix} a\\b\\c \end{matrix}\right) = \sum_{\textbf x\in S} \left(\begin{matrix}x\cdot I(\textbf x) \\ y\cdot I(\textbf x) \\I(\textbf x)\end{matrix}\right) \tag{7}

如果它不对劲（不满秩），说明至少有一个维度坍缩了，可以退一步尝试随便解两个变量，如果还坍缩，那就再退一步用常值代替插值方法。

结果

上图为2600个三角形的结果，可见该方法的显著缺陷是主次不分，本来就虚焦的背景却堆了大量的三角形，而重要的人像由于特殊摄影手法造成的低对比度而变得非常概括。

但毕竟只是个课程大作业，所以大概不会有后续跟进了。我们的代码在Github上开源。