写在前面
好久没有更新shadertoy系列了，我万万没想到有童鞋还惦记着它。。。之前说过希望可以一周更新一篇，现在看来是不怎么可能了，一个月更新一篇的希望比较大（不要再相信我了。。。）
我把之前实现的这个系列上传到了GitHub（https://github.com/candycat1992/Shadertoy_Lab）上，有兴趣的可以去下载下来。当然，也希望有网友可以一起贡献这个项目。
GitHub上这个项目大部分灵感来源于shadertoy（https://www.shadertoy.com），也有一些是配合博客里的一些文章讲解的，也有一些是在原shadertoy里面的例子扩展而来的。总之，每个lab我都会在README里面给出相关的参考资料。
项目链接：https://github.com/candycat1992/Shadertoy_Lab
好啦，我们来看一下这一篇里面要讲的例子。如同题目所讲，我们的目标是模拟水彩风格的效果。当然，这里实现的只是简化版本后的实现，我们只实现了渲染部分。
参考资料：
[1] https://www.shadertoy.com/view/XdSSWd
[2] Curtis C J, Anderson S E, Seims J E, et al. Computer-generated watercolor[C]// Proceedings of the 24th annual conference on Computer graphics and interactive techniques. ACM Press/Addison-Wesley Publishing Co., 1997.
论文研讨：Computer-Generated Watercolor
这个例子来源于一篇著名的论文，也就是1997年的Computer-Generated Watercolor。年代虽然很久远了，但是这篇论文开启了用计算机模拟水彩画的先河，后面陆陆续续又有很多论文被发表出来，但几乎都可以看到这篇论文的影子。
这篇论文主要可以分为四个部分：
- 首先，描述了水彩颜料的物理性质，并从艺术角度给出了一些水彩画的风格特性；
- 给出了如何模拟这些特性的方法；
- 具体描述了对水彩和颜料（pigment）的物理模拟算法；
- 描述了如何渲染这些颜料。
而本文其实只是实现了最后一个部分，在本节后面的内容里，我会简略介绍下论文里其他三个方面的内容。如果读者对这方面研究有兴趣的话，还是强烈建议去阅读原论文。
水彩的物理属性
水彩画（watercolor paint ，也被简称为watercolor）是一种比较常见的艺术风格。一幅水彩画涉及到了两种材质：
水彩纸（watercolor paper）。它并不是由木材制作而成的，而是通过把亚麻布或者棉花捣碎成细小的纤维的来的。这种材质非常容易吸收液体，为了防止颜料迅速蔓延，因此还给这些纸张进行上浆（sizing）。
颜料（pigment）。这是一种固体材质，由很多很小的单独的粒子组成。这些水彩颜料通常由0.05到0.5微米的粉末构成，它们可以渗透水彩纸，但一旦附着在纸上，扩散速度就会下降。
除此之外，水彩画有一些特点，例如：
干笔画（Dry brush）。如果使用较干的画笔画在粗糙的纸上，那么会出现一些不规则的空隙和粗糙的边界效果。
边界颜色较深（Edge darkening）。如果使用较湿的画笔画在较干的纸面上，在纸的浆料和水的表面张力的作用下，颜料不会继续扩散，并在边缘处留下一圈颜色更深的沉淀痕迹。
模拟
在论文中，作者提出使用三个图层来模拟水彩画中颜料的流动：
第一层是shallow-water layer。在这一层中，水和颜料会在纸张表面扩散流动。
第二层是pigment-deposition layer。在这一层中，颜料会沉淀进入和释放出纸张。
第三层是capillary layer。在这一层中，被纸张吸收的水会通过毛细管作用被继续扩散。（这一层仅仅用于模拟水彩画的回流效果。）
在模拟时，作者使用了很多参数来控制模拟效果，例如颜料的扩散速度、画笔压力、纸张的高度、颜料密度、液体饱和度、液体容量等等。
关于纸张的模拟，作者使用了一种简单的模型，即高度场的方法，并使用了Perlin噪声（Ken Perlin. An image synthesizer. In SIGGRAPH ’85 Proceedings, pages 287–296. July 1985.）和Worley的多孔纹理（Steven P. Worley. A cellular texturing basis function. In SIGGRAPH ’96 Proceedings, pages 291–294. 1996.）来生成。这种方法非常常见。
算法
有了上述这些参数之后，就可以进行算法模拟的部分。主循环部分在每个时间步内，会进行四个计算步骤：
在shallow-water layer移动液体（Move Water）。
在shallow-water layer移动颜料（Move Pigment）。
在pigment-deposition layer传递颜料（Transfer Pigment）。这一步会模拟颜料的吸收和释放。
在capillary layer模拟毛细流动（Simulate Capillary Flow）。这一步会模拟回流现象等。
具体的算法还是要参考论文，本文不涉及这些算法的实现。
渲染
以上的内容只是为了完整性，而与这篇博客相关的只有渲染部分。
当经过上面的算法后，我们可以得到每个区域的颜料厚度。
作者使用了Kubelka-Munk（KM）模型来渲染颜料。在论文中，作者为每个颜料指定了两个系数：吸收系数（absorption coefficients）K和散射系数（scattering coefficients）S。K和S都是三维属性，分别表示颜料吸收和散射的能量。
指定颜料的光学属性
虽然K和S系数通常是经验决定的，但作者允许让用户来指定：通过选择希望的“unit thickness”（单位厚度）的颜料在黑白背景下的外观来决定。具体方法是，给定用户选择的两个RGB颜色
作者在论文里给出了一些计算出来的不同样色、不同属性颜料的KS系数。如下图所示（图片来源《Computer-Generated Watercolor》）：
这些颜料是不同类型的，例如：
* 不透明颜料（Opaque paints）。类似Indian Red（上图中的b），在白色和黑色区域都有类似的颜色。这种颜料都具有高散射、高吸收的属性。
透明颜料（Transparent paints）。类似Quinacridone Rose（上图中的a），在白色背景下有颜色，在黑色背景下几乎是黑色的。这种颜料的scattering波长都很低，而absorption分量很高，并和它们的颜色是补集。
干涉颜料（Interference paints）。类似Interference Lilac（上图中的l），在白色背景下几乎是白色的，而在黑色背景下是有颜色的。
光学的颜料层混合
一旦给定了一个一定厚度x的颜料层以及它的散射和吸收系数S和K，我们就可以按下面的公式计算该颜料层的反射比R和透射比T。
对于两个相邻的层，我们可以按下面公式来计算合成后的颜料层的R和T：
Shader的实现
下面的内容会解释如何使用Unity Shader来实现上面的渲染部分。
从上面的渲染算法中可以看出，实际上渲染部分只涉及到了每个区域的颜料厚度x以及颜料的系数K和S。在下面的实现中，我们使用了论文中提供的一系列K和S。为此，我们在shader中定义了如下变量：
// Table of pigments
// from Computer-Generated Watercolor. Cassidy et al.
// K is absorption. S is scattering
// a
#define K_QuinacridoneRose vec3(0.22, 1.47, 0.57)
#define S_QuinacridoneRose vec3(0.05, 0.003, 0.03)
// b
#define K_IndianRed vec3(0.46, 1.07, 1.50)
#define S_IndianRed vec3(1.28, 0.38, 0.21)
// c
#define K_CadmiumYellow vec3(0.10, 0.36, 3.45)
#define S_CadmiumYellow vec3(0.97, 0.65, 0.007)
// d
#define K_HookersGreen vec3(1.62, 0.61, 1.64)
#define S_HookersGreen vec3(0.01, 0.012, 0.003)
// e
#define K_CeruleanBlue vec3(1.52, 0.32, 0.25)
#define S_CeruleanBlue vec3(0.06, 0.26, 0.40)
// f
#define K_BurntUmber vec3(0.74, 1.54, 2.10)
#define S_BurntUmber vec3(0.09, 0.09, 0.004)
// g
#define K_CadmiumRed vec3(0.14, 1.08, 1.68)
#define S_CadmiumRed vec3(0.77, 0.015, 0.018)
// h
#define K_BrilliantOrange vec3(0.13, 0.81, 3.45)
#define S_BrilliantOrange vec3(0.009, 0.007, 0.01)
// i
#define K_HansaYellow vec3(0.06, 0.21, 1.78)
#define S_HansaYellow vec3(0.50, 0.88, 0.009)
// j
#define K_PhthaloGreen vec3(1.55, 0.47, 0.63)
#define S_PhthaloGreen vec3(0.01, 0.05, 0.035)
// k
#define K_FrenchUltramarine vec3(0.86, 0.86, 0.06)
#define S_FrenchUltramarine vec3(0.005, 0.005, 0.09)
// l
#define K_InterferenceLilac vec3(0.08, 0.11, 0.07)
#define S_InterferenceLilac vec3(1.25, 0.42, 1.43)对于颜料厚度，则是基于distance field的方法，再通过一些计算来模拟Edge darkening的效果。
我们首先来看颜料层的反射比R和透射比T。代码如下：
// Kubelka-Munk reflectance and transmitance model
void KM(vec3 K, vec3 S, float x, out vec3 R, out vec3 T) {
vec3 a = (K + S) / S;
vec3 b = sqrt(a * a - vec3(1.0));
vec3 bSx = b * S * vec3(x);
vec3 sinh_bSx = my_sinh(bSx);
vec3 c = a * sinh_bSx + b * my_cosh(bSx);
R = sinh_bSx / c;
T = b / c;
}它的输入有三个：该区域颜料的吸收系数K和散射系数S，颜料厚度x。输出该区域的反射比R和透射比T。有了上一节的公式，上面的代码就很简单了，就是带公式而已。
另一个公式是用于混合两个颜料层。相关代码如下：
// Kubelka-Munk model for optical compositing of layers
void CompositeLayers(vec3 R0, vec3 T0, vec3 R1, vec3 T1, out vec3 R, out vec3 T) {
vec3 tmp = vec3(1.0) / (vec3(1.0) - R0 * R1);
R = R0 + T0 * T0 * R1 * tmp;
T = T0 * T1 * tmp;
}它的输入是两个颜料层的反射比和透射比，输出合成层的反射比和透射比。同样，上面的代码也是带公式而已。
至此，我们想要渲染出来画面还需要提供一个参数，就是KM函数中的颜料厚度x。在论文中，这个颜料厚度是通过一系列算法计算得到的。但在我们的实现中，我们简化了这一步，而使用基于distance field的方法来计算厚度。要模拟一定效果的水彩风格，我们需要模拟它的一些特性，例如Dry-brush和Edge darkening。我们这里选择一种取巧的方法，利用了噪声（模拟粗糙的边界效果）和一些数学计算（来模拟Edge darkening效果），而非原文中复杂的算法。
在我们的实现中，我们是在fragment shader中渲染图形的，也就是，我们处理的单位是逐像素的。当渲染一个图形时，我们需要几个步骤：
给定渲染区域的位置pos。为了模拟水彩画粗糙的边缘效果，我们会使用噪声函数来对屏幕坐标进行一些的噪声处理。例如：
vec2 uv = fragCoord.xy / iResolution.xy ... pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.02 * noise2d(uv * vec2(0.1)));uv是对屏幕坐标处理后得到的xy范围都在(0, 1)的坐标。为了计算当前位置的坐标，我们首先针对屏幕长宽进行处理，使得得到的坐标在x方向上的范围是(0, 1)，y方向上的范围是(0, height/width)。然后，对结果添加了噪声处理。
噪声函数noise2d的代码如下：
// Simple 2d noise fbm (Fractional Brownian Motion) with 3 octaves float Noise2d(vec2 p) { float t = texture2D(iChannel0, p).x; t += 0.5 * texture2D(iChannel0, p * 2.0).x; t += 0.25 * texture2D(iChannel0, p * 4.0).x; return t / 1.75; }这是一种非常简单的噪声实现。主要通过对一张噪声纹理采样，并使用了三层的octaves。这些内容可以在Perlin噪声的相关内容中找到，例如这里（http://freespace.virgin.net/hugo.elias/models/m_perlin.htm）和这里（http://stackoverflow.com/questions/16999520/in-need-of-fractional-brownian-noise-fbm-vs-perlin-noise-clarification）。我们简单解释一下，不同的octave表示不同的频率和振幅噪声，通过组合不同频率和振幅的噪声，我们就可以得到一个Perlin噪声。
在得到了区域坐标pos后，我们需要把它代入distance field的计算，得到距离值dist。例如：
dist = DistanceCircle(pos, vec2(0.2, 0.55), 0.08);DistanceCircle函数代码如下：
float DistanceCircle(vec2 pos, vec2 center, float radius) { return 1.0 - distance(pos, center) / radius; }它会计算pos相对于圆心在center、半径为radius的圆的距离，返回值 > 0时表示在圆内，返回值 < 0时表示在圆外。
类似的距离计算函数还有DistanceLine（对应画直线）、DistanceSegment（对应画线段）和DistanceMountain（对应画一座由正弦函数得到的山脉）等。
在得到了距离值后，我们由此来判断一个点是否需要绘制水彩。但为了后面的渲染，我们还需要把这个距离值转换成颜料厚度，这是通过BrushEffect函数得到的。例如：
float circle = BrushEffect(dist, 0.2, 0.1);BrushEffect函数代码如下：
// Simulate edge darkening effect // Input: dist < 0 outer area, dist > 0 inner area float BrushEffect(float dist, float x_avg, float x_var) { // Only when abs(dist) < 1.0/10.0, x > 0.0 // Means that the edges have more thickness of pigments float x = max(0.0, 1.0 - 10.0 * abs(dist)); x *= x; x *= x; return (x_avg + x_var * x) * smoothstep(-0.01, 0.002, dist); }BrushEffect不仅会把距离值dist变换到颜料厚度，也会负责模拟Edge darkening效果。它的输入是上一步计算而得的dist（dist < 0表示在渲染图形的外部，dist > 0表示在内部），以及平均颜料厚度x_avg和边缘厚度变化x_var。计算过程是：
第一行首先根据dist计算初始的边缘颜料厚度x，它的范围是(0, 1)。当dist的绝对值小于1/10时（即靠近边界处），x大于0；否则x等于0。我们也可以调整公式中的参数10，值越小，Edge darkening的范围就越广。
后面两行进一步处理边缘颜料厚度x，这是通过自乘两次得到的。这样可以进一步快读收紧Edge darkening的范围。
计算返回值即颜料的厚度。首先我们通过smoothstep函数来控制厚度的整体变化，具体是，当dist小于-0.01时，返回0，当大于0.002返回1，否则返回0到1之间的值。-0.01和0.002的选择并不是完全任意的，我们一般选择一正一负来处理边界，同时正数（这里是0.002）的数值通常要小于负数的绝对值（这里是|-0.01|），这是为了让来模拟出颜料在边界处的扩散速度非线性下降的效果。然后，我们把该值和(x_avg + x_var * x)的结果相乘。其中，x_avg是渲染图形内部绝大多数区域的颜料厚度，而x_var用于控制边界处的颜料厚度（因为边界处的颜料厚度要大于内部），x_var越大，边界出的Edge darkening效果越明显。在我们的实现中，一般取x_avg为0.2，取x_var为0.1。当然，如果我们想要模拟出粗糙感，也可以传噪声进去，例如：
float mountains = BrushEffect(dist, 0.2, 0.3 * Noise2d(uv * vec2(0.1)));至此，我们就得到了该区域的颜料厚度。同样，为了简单的模拟颜料不均匀分布的特点，我们也可以进一步对结果值进行噪声处理。例如：
mountains *= 0.65 + 0.35 * Noise2d(uv * vec2(0.2));注意其中的系数0.65和0.35，它们的和需要是1。如果把0.35的调大、把0.65调小，粗糙感就越强烈。
最后，我们只需要把颜料厚度，和选择的KS系数传递给KM函数得到该颜料层的反射比和透射比。如果需要和之前的颜料层混合，再代入CompositeLayers函数混合即可。例如：
KM(K_HansaYellow, S_HansaYellow, circle, R1, T1); CompositeLayers(R0, T0, R1, T1, R0, T0);
实现效果
shadertoy中原作者的绘制结果我在Unity中重现后是下面这样的（调整了一些参数）：
上述场景的绘制代码如下：
///
/// First Scene
///
// Background
float background = 0.1 + 0.1 * Noise2d(uv * vec2(1.0));
KM(K_CeruleanBlue, S_CeruleanBlue, background, R0, T0);
pos = uv + vec2(0.04 * Noise2d(uv * vec2(0.1)));
dist = DistanceMountain(pos, 0.5);
float mountains = BrushEffect(dist, 0.2, 0.3 * Noise2d(uv * vec2(0.1)));
mountains *= 0.45 + 0.55 * Noise2d(uv * vec2(0.2));
KM(K_HookersGreen, S_HookersGreen, mountains, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.02 * Noise2d(uv * vec2(0.1)));
dist = DistanceCircle(pos, vec2(0.2, 0.55), 0.08);
float circle = BrushEffect(dist, 0.2, 0.2);
KM(K_HansaYellow, S_HansaYellow, circle, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);我在原shader的基础上做了一些扩展，例如给出了原论文中所有样例的KS系数，给出了更过的距离计算函数。通过这些的组合，可以得到更多的效果。例如：
上述场景的绘制代码如下：
///
/// Second Scene
///
// Background
float background = 0.1 + 0.2 * Noise2d(uv * vec2(1.0));
KM(K_HansaYellow, S_HansaYellow, background, R0, T0);
// Edge roughness: 0.04
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.04 * Noise2d(uv * vec2(0.1)));
dist = DistanceCircle(pos, vec2(0.5, 0.5), 0.15);
// Average thickness: 0.2, edge varing thickness: 0.2
float circle = BrushEffect(dist, 0.2, 0.2);
// Granulation: 0.85
circle *= 0.15 + 0.85 * Noise2d(uv * vec2(0.2));
KM(K_CadmiumRed, S_CadmiumRed, circle, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);
// Edge roughness: 0.03
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.03 * Noise2d(uv * vec2(0.1)));
dist = DistanceCircle(pos, vec2(0.4, 0.3), 0.15);
// Average thickness: 0.3, edge varing thickness: 0.1
circle = BrushEffect(dist, 0.3, 0.1);
// Granulation: 0.65
circle *= 0.35 + 0.65 * Noise2d(uv * vec2(0.2));
KM(K_HookersGreen, S_HookersGreen, circle, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);
// Edge roughness: 0.02
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.02 * Noise2d(uv * vec2(0.1)));
dist = DistanceCircle(pos, vec2(0.6, 0.3), 0.15);
// Average thickness: 0.3, edge varing thickness: 0.2
circle = BrushEffect(dist, 0.3, 0.2);
// Granulation: 0.45
circle *= 0.55 + 0.45 * Noise2d(uv * vec2(0.2));
KM(K_FrenchUltramarine, S_FrenchUltramarine, circle, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);
// Opaque paints, e.g. Indian Red
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.02 * Noise2d(uv * vec2(0.3)));
dist = DistanceSegment(pos, vec2(0.2, 0.1), vec2(0.4, 0.25), 0.03);
float line = BrushEffect(dist, 0.2, 0.1);
KM(K_IndianRed, S_IndianRed, line, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);
// Transparent paints, e.g. Quinacridone Rose
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.02 * Noise2d(uv * vec2(0.2)));
dist = DistanceSegment(pos, vec2(0.2, 0.5), vec2(0.4, 0.55), 0.03);
line = BrushEffect(dist, 0.2, 0.1);
KM(K_QuinacridoneRose, S_QuinacridoneRose, line, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);
// Interference paints, e.g. Interference Lilac
pos = uv * vec2(1.0, iResolution.y / iResolution.x) + vec2(0.02 * Noise2d(uv * vec2(0.1)));
dist = DistanceSegment(pos, vec2(0.6, 0.55), vec2(0.8, 0.4), 0.03);
line = BrushEffect(dist, 0.2, 0.1);
KM(K_InterferenceLilac, S_InterferenceLilac, line, R1, T1);
CompositeLayers(R0, T0, R1, T1, R0, T0);注意到上面对参数的调整和不同图形效果的区别，例如边界颜色更深、边界粗糙感和整体颗粒感等等。
完整的代码读者可以在https://github.com/candycat1992/Shadertoy_Lab中的WaterColorScene找到。
写在最后
本文实现了水彩风格的渲染部分。对于颜料厚度的计算则是通过简单的数学计算来模拟的，当然，这样得到的效果也并不真实。想要得到更加真实的效果，需要配合更复杂的算法，具体可以参见上面的论文及其它发表的论文。
这篇文章只是抛砖引玉，从它的实现我们可以学到KM模型在实时渲染中的实现，以及噪声的简单应用。在本文中，我们的实现都是基于distance field的方法，也就是说我们画出的图形其实都是用数学表达式计算而得的。读者可以添加更过的函数来绘制更复杂的图形。如果想要实现那种用户交互的应用，也可以使用其它方法来计算颜料厚度。
最后，希望这篇文章可以对大家有所帮助~