Creation by Silexars
Creation by Silexars

前言

最近无聊闲逛ShaderToy时，无意间发现Danguafer大佬的一个作品，上面是这个作品的截图。
五彩的耀斑从中心四散而开，非常的amazing！

这时我心想，肯定是大佬用复杂的算法实现的，学不来学不来。
当我刚要按下Ctrl+W时，我瞄了一眼代码，发现这个复杂的图形竟然只是不到30行代码实现的！
这让我更amazing了，这个算法究竟是使用了什么原理呢？
我决定一探究竟。

废话

下面文章内容将围绕这30行代码展开，阅读大概需要10分钟。
要看懂这篇文章，首先你需要一(亿)点点线代基础和图形学知识。
哈哈，不玩梗了，其实也不用太多基础知识。
考虑到阅读体验，我的文章使用直观的可视化图形讲解，尽量不会出现大片的数学公式，放心食用。

文献和声明

ShaderToy原文地址：https://www.shadertoy.com/view/XsXXDn
Danilo Guanabara：http://www.pouet.net/prod.php?which=57245
GLSL文档推荐：https://github.com/wshxbqq/GLSL-Card

本文非商业文章，且为了方便大家阅读，我将源代码复制到下面，供大家参考。
其著作权归原作者Danguafer所有，若想商用请联系原作者。
按照大佬在原文1，2行的规定，我上面放出Danguafer作者文章地址。

本文全部动画皆为原创内容，若要转载请联系我，并著名出处。
作者不是专业研究数学的，如有错误请各位谅解。

本文作者联系方式：mrkbear@qq.com

Danguafer大佬的源代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// http://www.pouet.net/prod.php?which=57245
// If you intend to reuse this shader, please add credits to 'Danilo Guanabara'

#define t iTime
#define r iResolution.xy

void mainImage( out vec4 fragColor, in vec2 fragCoord ){
	vec3 c;
	float l,z=t;
	for(int i=0;i<3;i++) {
		vec2 uv,p=fragCoord.xy/r;
		uv=p;
		p-=.5;
		p.x*=r.x/r.y;
		z+=.07;
		l=length(p);
		uv+=p/l*(sin(z)+1.)*abs(sin(l*9.-z*2.));
		c[i]=.01/length(abs(mod(uv,1.)-.5));
	}
	fragColor=vec4(c/l,t);
}

正式开始

首先以我们的最大能力去试着读一下这段代码。

从4，5行可以得知，这个Shader程序的全局变量：

1
2
#define t iTime            //使用t变量接收当前时间，随着程序执行，时间不断增大 
#define r iResolution.x    //使用向量r接收，画布的尺寸数据

从7行可以得知，这个Shader程序的输入变量：

1
2
3
4
void mainImage(			//Shader的主函数，相当于C语言的Main
	out vec4 fragColor,	//输出四维向量 片段颜色
	in vec2 fragCoord	//输入二维向量 片段位置
)

简单分析一下程序的逻辑结构：

定义vec3 c作为输出结果
定义float l存储长度，float z备份时间

接下来的for循环遍历c的每一个(色彩)维度，以0.7为相位计算，每个维度颜色的值。
最后将c除上l，在与时间拼合，作为输出值。
因为每个维度在time上差0.7所以看上去是彩色的。

简化程序

现在我们已经对整个程序大体上有一定了解了。

想要深入了解其中的核心算法，就要一层一层简化这个程序，就像剥洋葱一样，把核心算法提取出来，最后你会感动都落泪。
这个程序要遍历三个颜色维度，而咱们连一个都玩不明白，别说三个了。

下面我们将其一个颜色维度的核心算法提取出来，得到这样的程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// 将坐标空间缩放到 ([-0.5 ~ 0.5], [-0.5 ~ 0.5])
vec2 pos = v_texcoord - .5;

// 保持 Y 坐标不变，将 X 坐标空间缩放到画布比例
pos.x *= resolution.x / resolution.y;

// 计算向量 片段位置距离屏幕中心的距离
float l = length(pos);

// 看不懂的神仙代码1
pos = v_texcoord + pos / l * (sin(time) + 1.) * abs(sin(l * 9. - time * 2.));

// 看不懂的神仙代码2
float res = .01 / length(abs(mod(pos, 1.) -.5));

// 返回颜色向量
fragColor = vec4(vec3(res / l), 1);

我们把上面的代码放到 HLSL 程序中跑一遍，看看结果。

低配Creation by Silexars
低配Creation by Silexars

果然，我们得到了一个低配版的Creation by Silexars！
它的效果和原程序一模一样，但是它只有一个色彩维度。
这样，我们把原程序13行的核心算法简化为6行！！！

什么，你还是觉得它很复杂？
那咱们再删掉一行。

我们来看第3行：

1
pos.x *= resolution.x / resolution.y;

假如让resolution.x与resolution.y相等，我们就能删掉这行了。
这样做同时也会降低之后建模的复杂度。
让resolution.x与resolution.y相等也就是让图形绘制在正方形画布上。

正方形Creation by Silexars
正方形Creation by Silexars

来，安排。

抽象成数学公式

数学公式
数学公式

为了方便我们进一步骤研究，我们将之前看不懂的两行神仙代码，抽象成简单的数学公式。
可以看到是结果是分4步计算的：

  1. 计算当前片段与远点之间的距离。
  2. 计算两个奇怪的正弦波的乘积，其中一个正弦波经过绝对值处理。
  3. 将当前向量转换为标准向量并与第2步的值相乘，在加上当前向量…
  4. 不废话了，算法在公式里

使用matlab建模

为了直观的看到，前面公式中每一步的函数波形变化，我在matlab中实现了这个算法。
并尝试用matlab绘制出和GPU一样的图形。

经过不断的尝试，使用matlab实现的算法如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
% 时间
time = 0;
% 网格精细度
sp = 0.003;

[X,Y] = meshgrid(-0.5:sp:0.5, -0.5:sp:0.5);

L = (X.^2 + Y.^2).^(1/2);

% 这里我们称之为“第一阶段”
Xb = X + .5 + X ./ L .* (sin(time) + 1.) .* abs(sin(L * 9. - time * 2.));
Yb = Y + .5 + Y ./ L .* (sin(time) + 1.) .* abs(sin(L * 9. - time * 2.));

% 这里我们称之为“第二阶段”
Xb = abs(mod(Xb, 1.) -.5);
Yb = abs(mod(Yb, 1.) -.5);

NL = .01 ./ (Xb.^2 + Yb.^2).^(1/2) ./ L;

surf(X, Y, NL);

注意这里我们仅仅绘制time=0时的图像。

失败品
失败品

结果令我们非常失望，一片汪洋大海，中间一柱擎天。
难道是算法出了问题嘛？
后来，经过我仔细观察才发现，原来是中间靠近(x=0,y=0)的点的数值太高了，导致我们无法看清其他位置的细节。
而实际上我们的图像已经成功绘制了。

为了解决这个问题，我们在绘制图形前，将结果进行开8次根号处理，这样原本非常高的数值相对变得小了很多。
最好还是要修改坐标轴刻度排列方式为非线性，但是不知道为什么，这样做在绘制动画时会出现一些小问题，也许今后有机会我会修复它。

1
NL = NL .^ (1/8);

处理后绘制的图像：

成功品
成功品

这3张图片里，左上角、右上角、下面的图分别时是在time=0时GPU生成的图像、matlab绘制的图像、matlab绘制的三维曲面图。

这里我们可以看到，我们的matlab模型，和GPU图像完全拟合了。
到这里我们成功了50%了，而剩下的50%，就是使用matlab模型一步一步研究算法各个时期的波形图。

分析L的图像

1
L = (X.^2 + Y.^2).^(1/2);

不难看出，L与time是无关的，且L代表当前片段坐标与原点的距离。
我们可以很轻松的在脑海中想象出L关于X,Y的图像，大概像一个大漏斗。

什么？你想象不出来？
哎哎哎，客官别走，我这就给您画一下。

漏斗
漏斗

分析两段正弦波图像

下面我们来分析，整个算法中，两段正弦波，以及它们相乘的图像：

1
2
3
wave1 = sin(time) + 1.;
wave2 = abs(sin(L * 9. - time * 2.));
wave3 = wave1 .* wave2;

首先wave1的图像长这样，不多解释了：

wave1
wave1

然而wave2的图像就没有那么简单了。
首先wave2与L相关，同时又与time相关。

首先来看L*9是把“漏斗”沿着Z轴放大9倍。
而-time * 2的含义是：“漏斗”不断的沿着Z轴负方向移动。
如果你看到上面两句话，并在脑海里想象出一个不断下坠的漏斗的图像，那么恭喜你，你理解了。

下面我们来看wave2图像。
差不多将刚刚想象的“不断下坠的漏斗”融入“波的图像”，并随着time的变化不断波动的样子。
这里要注意，因为绝对值的影响，z>0的部分要被翻折到z=0平面上方。
当然，能想象出这个，并不容易。

为此，我在matlab中使用timer制作了一个小动画帮助你想象这个画面。

wave2
wave2

如果你想象出的画面大概与这个吻合，我不得不承认你的空间想象能力比我厉害。

最后我们来想象一下wave3的图像，将上面wave2和wave1的图像相乘。
这个还是比较好理解的，差不多是，上面波的高度随着wave1不断的从0到2从2到0。

来吧，老样子，上图：

wave3
wave3

注：wave2的gif图2.7M大小，wave3的gif图8.3M大小，注意流量！！！

分析X./L与Y./L的图像

X./L与Y./L的区别仅仅在于数据的方向不同，所以他们图像应该非常相似。
所以我们抛下一个稍后研究，下面我们来单独分析X./L的图像。

别看Y./L看似简单，实际上线性函数和双曲线的复合在乘上一条线性函数，结果绝对出乎你的意料。
我们将Y./L拆解成2部分，分别研究。

1
2
dp1 = L .^ (-1);
dp2 = X .* dp1;

首先来看dp1的图像，在看结果之前，我们不妨在脑海中想象一下这个图像是什么样的。
首先L的图像是个漏斗，代表它的值是从(0,0)点到任意位置的过程是线性增长的。
是不是很像函数y=x在x>0时的图像？
由此我们能想到y=x .^ (-1)的图像就是我们熟悉的双曲线y=1/x。
我们大胆的猜测dp1的图像应该是，y=1/x在x>0时的图像，绕着Z轴旋转得到的图形。

看完上面的几句话，并在脑海里想象出一个“锥子”的图像，那么恭喜你，你理解了。

dp1
dp1

dp2的图像：

dp2
dp2

第一阶段的图像

1
Xb = X + .5 + X ./ L .* (sin(time) + 1.) .* abs(sin(L * 9. - time * 2.));

将上面的图像，按照公式运算后，得出第一阶段的图像：

dp2
dp2

第二阶段的图像

1
Xb = abs(mod(Xb, 1.) -.5);

通过上面的公式可以看出，二阶段做了3件事：

  1. 对于1取模
  2. -0.5
  3. 取绝对值

对于1取模会把不在0到1之间值映射到0到1之间。
这样做可以限制图像值域到(0,1)之间。

减0.5后取绝对值，会将小于0.5的部分，翻折到z>0的空间中。

最终得到的第二阶段图像：

第二阶段图像
第二阶段图像

最后一步图像

1
2
3
dis = (Xb.^2 + Yb.^2).^(1/2);
dos = dis ^ (-1);
NL = .01 * dos ./ L;

终于到最后一步了，确认存活？

dis在计算第二阶段产生向量的长度。
dos是y=1/x与dis复合产生的。

最终图像

经过一系列运算，终于…

这是最终得到的图像，它仿佛像一个艺术作品：

最终图像
最终图像

注意：最后这个模型，切除了z>1的部分，因为fragColor接收的有效值在0到1之间。

写在最后

数学也是一种艺术，它同样有鲜艳的色彩，同样有美妙的声音，同样有动感的画面。

非常感谢你能看到这里！！