【ShaderToy】水彩画

写在前面

好久没有更新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。

年代尽管非常久远了，可是这篇论文开启了用计算机模拟水彩画的先河，后面陆陆续续又有非常多论文被发表出来。但差点儿都能够看到这篇论文的影子。

这篇论文主要能够分为四个部分：

1. 首先，描写叙述了水彩颜料的物理性质。并从艺术角度给出了一些水彩画的风格特性；
2. 给出了怎样模拟这些特性的方法；
3. 详细描写叙述了对水彩和颜料（pigment）的物理模拟算法；
4. 描写叙述了怎样渲染这些颜料。

而本文事实上仅仅是实现了最后一个部分。在本节后面的内容里，我会简略介绍下论文里其他三个方面的内容。假设读者对这方面研究有兴趣的话，还是强烈建议去阅读原论文。

水彩的物理属性

水彩画（watercolor paint ，也被简称为watercolor）是一种比較常见的艺术风格。一幅水彩画涉及到了两种材质：

• 水彩纸（watercolor paper）。它并非由木材制作而成的，而是通过把亚麻布或者棉花捣碎成细小的纤维的来的。这样的材质非常easy吸收液体，为了防止颜料迅速蔓延，因此还给这些纸张进行上浆（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.）来生成。这样的方法非经常见。

算法

有了上述这些參数之后，就能够进行算法模拟的部分。主循环部分在每一个时间步内。会进行四个计算步骤：

1. 在shallow-water layer移动液体（Move Water）。

2. 在shallow-water layer移动颜料（Move Pigment）。

3. 在pigment-deposition layer传递颜料（Transfer Pigment）。

   这一步会模拟颜料的吸收和释放。

4. 在capillary layer模拟毛细流动（Simulate Capillary Flow）。这一步会模拟回流现象等。

详细的算法还是要參考论文，本文不涉及这些算法的实现。

渲染

以上的内容仅仅是为了完整性，而与这篇博客相关的仅仅有渲染部分。

当经过上面的算法后，我们能够得到每一个区域的颜料厚度。

作者使用了Kubelka-Munk（KM）模型来渲染颜料。

在论文中。作者为每一个颜料指定了两个系数：吸收系数（absorption coefficients）K和散射系数（scattering coefficients）S。K和S都是三维属性，分别表示颜料吸收和散射的能量。

指定颜料的光学属性

尽管K和S系数一般是经验决定的，但作者同意让用户来指定：通过选择希望的“unit thickness”（单位厚度）的颜料在黑白背景下的外观来决定。详细方法是，给定用户选择的两个RGB颜色Rw（在白色背景下的颜色）和Rb（在黑色背景下的颜色），K和S系数能够靠以下的等式来得到：

S=1b⋅coth−1(b2−(a−Rw)(a−1)b(1−Rw))K=S(a−1)其中。a=12(Rw+Rb−Rw+1Rb)b=a2−1‾‾‾‾‾‾√

作者在论文里给出了一些计算出来的不相同色、不同属性颜料的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=sinhbSxcT=bc，其中c=asinhbSx+bcoshbSx

对于两个相邻的层，我们能够按以下公式来计算合成后的颜料层的R和T：

R=R1+T21R21−R1R2T=T1T21−R1R2

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。有了上一节的公式，上面的代码就非常easy了，就是带公式而已。

还有一个公式是用于混合两个颜料层。

相关代码例如以下：

    // 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中渲染图形的。也就是，我们处理的单位是逐像素的。当渲染一个图形时，我们须要几个步骤：

1. 给定渲染区域的位置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;
   }

   这是一种非常easy的噪声实现。

   主要通过对一张噪声纹理採样。并使用了三层的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噪声。

2. 在得到了区域坐标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（相应画一座由正弦函数得到的山脉）等。

3. 在得到了距离值后。我们由此来推断一个点是否须要绘制水彩。但为了后面的渲染，我们还须要把这个距离值转换成颜料厚度，这是通过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。

   计算过程是：

   1. 第一行首先依据dist计算初始的边缘颜料厚度x，它的范围是(0, 1)。

      当dist的绝对值小于1/10时（即靠近边界处）。x大于0。否则x等于0。我们也能够调整公式中的參数10，值越小，Edge darkening的范围就越广。

   2. 后面两行进一步处理边缘颜料厚度x。这是通过自乘两次得到的。这样能够进一步快读收紧Edge darkening的范围。

   3. 计算返回值即颜料的厚度。首先我们通过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调小，粗糙感就越强烈。

4. 最后。我们仅仅须要把颜料厚度，和选择的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的方法，也就是说我们画出的图形事实上都是用数学表达式计算而得的。

读者能够加入更过的函数来绘制更复杂的图形。假设想要实现那种用户交互的应用。也能够使用其他方法来计算颜料厚度。

最后，希望这篇文章能够对大家有所帮助~