IFS 分形揭秘

2024-06-16

IFS 方法
轨道着色
着色器编程基础
空间折叠
实战演示

本文整理自我 2024/06/14 在上海科技大学数学所的一个小报告，标题是「GPU 涂鸦与数学可视化」。这个报告是我和上科大数学所的陈浩老师、Abdelaziz Nait Merzouk 合作的双曲反射群画展期间的一个助兴小节目（感谢陈浩老师一己之力将画展办起来，并邀请我去上海逛了一圈）。我保留了报告的技术内容，略去了关于 demoscene 和分形文化的部分。

在 Shadertoy 上有很多效果酷炫，但是代码非常短的分形作品。我挑选了其中三个优秀的例子展示如下：

Ethereal by Kali	Apollonian fractal by Xor	Radiosity by Xor

然而代码短可不代表它们容易看懂。尤其是很多作者还喜欢故弄玄虚，把代码作了混淆处理以增加神秘感。对我来说，这种被人秀了一脸却没搞明白对面是怎么装的逼的感觉让人很不爽。当然我不是在抱怨，这种炫技的行为本身就是黑客文化的一部分，可以理解。后来分形玩的多了，我也慢慢明白了其中的奥妙，这次上科大之行是一次很好的机会，促使我把这些理解完整的写下来。

在这篇文章中，我将为大家揭示这些作品背后的奥秘。这些分形作品别看场景千变万化，其实都是用同一个套路制作出来的。这个套路可以简述为三步：

首先将像素的 2D 坐标映射为空间中的某个 3D 点 p；
然后用一个 fold 函数，即所谓的迭代函数系统（iterated function system，简称 IFS）反复作用在 p上，将 p 变换到空间中另一个位置。每次迭代结束后，生成一个颜色并将其添加到当前的颜色 color 上。
当达到一定的迭代次数后，终止迭代，返回 color 的值作为像素最终的颜色。

下面是具体步骤的解释。

IFS 方法

压缩映射.

设 \(f:\mathbb{R}^n\to\mathbb{R}^n\) 是一个映射，如果存在 \(0<r<1\) 使得对任何 \(x,y\in\mathbb{R}^n\) 有 \[d(f(x),f(y))\leq r\cdot d(x,y)\] 成立，我们就称 \(f\) 是一个压缩映射，\(r\) 是压缩比例。这里 \(d\) 是通常的 Euclidean 距离。

记 \(K(\mathbb{R}^n)\) 是 \(\mathbb{R}^n\) 中所有紧集组成的集合，可以证明 \(K(\mathbb{R}^n)\) 在 Hausdorff 度量下构成一个完备度量空间。读者不必关心 Hausdoff 度量的具体细节，只要知道它是一个度量，可以衡量两个紧集的接近程度。

设 \(f_1,\ldots,f_N\) 是 \(N\) 个压缩映射，\(f_i\) 的压缩比例是 \(r_i\)。定义映射 \(F:K(\mathbb{R}^n)\to K(\mathbb{R}^n)\) 如下： \[F(X) = f_1(X)\cup f_2(X)\cup\cdots\cup f_N(X),\quad X\in K(\mathbb{R}^n).\] 即 \(F\) 把 \(X\) 变成 \(N\) 个更小的集合。

可以证明 \(F\) 是空间 \(K(\mathbb{R}^n)\) 上的压缩映射，其压缩比例 \(r=\max\{r_1,\ldots,r_N\}\)。于是根据 Banach 不动点定理，存在唯一的紧集 \(A\subset\mathbb{R}^n\) 使得 \(A\) 是 \(F\) 的不动点： \[F(A) = A.\] 紧集 \(A\) 叫做 \(F\) 的极限集。\(A\) 是一个分形，它具有自相似的特征。

不仅如此，对 \(K(\mathbb{R}^n)\) 中的任何一点 \(K\)（\(K\) 是紧集），都有 \[\lim_{n\to\infty} F^n(K) =A.\]

我们以著名的 Sierpiński 三角形为例来说明这个过程。我们选择的三个压缩映射分别是

\[ \begin{align} f_1(x,y) &= (x/2, y/2),\\ f_2(x,y) &= (x/2, y/2) + (0, 1/2),\\ f_3(x,y) &= (x/2, y/2) + (1/2, 0).\\ \end{align} \]

这三个压缩映射的压缩比都是 1/2。初始的紧集 \(K_0\) 可以随便选，比如就取为一个圆：

\[K_0 = \img{/images/ifs/dot0.svg}{-1.75em}{}{4em}.\]

在第 1 次迭代后，它变成

\[K_1 = f_1(K_0)\cup f_2(K_0)\cup f_3(K_0) = \img{/images/ifs/dot1.svg}{-1.75em}{}{4em}.\]

在第 2 次迭代后，结果是

\[K_2 = f_1(K_1)\cup f_2(K_1)\cup f_3(K_1) = \img{/images/ifs/dot2.svg}{-1.75em}{}{4em}.\]

第 3 次迭代：

\[K_3 = f_1(K_2)\cup f_2(K_2)\cup f_3(K_2) = \img{/images/ifs/dot3.svg}{-1.75em}{}{4em}.\]

当迭代次数趋于无穷，就得到了 Sierpiński 分形：

\[\lim_{n\to\infty} K_n = A = f_1(A)\cup f_2(A)\cup f_3(A) = \img{/images/ifs/dot6.svg}{-1.75em}{}{4em}.\]

你可以很容易看出来为什么初始紧集的选择是不重要的：因为在压缩的过程中，任何紧集都会逐渐缩小到一个单点，所以 \(K_0\) 甚至取成一个点也是可以的。

压缩映射是无穷无尽的，所以 IFS 分形也是无穷无尽的。为了避免选择困难，我们一般只使用仿射和球反演这两种变换，通过组合它们来实现空间压缩。

注：这里的球反演变换是指将单位球的外部反演到内部，单位球内部保持不动的变换。虽然在单位球内部它并不压缩距离，但是只要和其它变换适当结合（必须有缩放），使得最终的复合变换是压缩的，就仍然可以得到分形。

然而，同时选择 \(N\) 个不同的变换 \(f_1,\ldots,f_N\)，还要让它们互相协调合作，共同生成漂亮的图案，还是太难做到了。我们后面会介绍，在着色器中绘制 IFS 是通过 \(F\) 的逆映射 \(F^{-1}\) 来实现的，\(F^{-1}\) 就是本文开头提到的 fold 函数，它是一种「空间折叠」操作，可以通过复合若干“折叠”函数来实现。就像古代的炼丹师会通过反复调整原料配方来寻找效果最佳的丹药一样，设计分形也可以通过调整 \(F^{-1}\) 中的折叠函数，并观察屏幕上显示的效果来实现。所以你根本无需关心 \(f_i\) 是什么！

轨道着色

我们希望给分形染上漂亮的颜色，这个染色应该满足如下的条件：

在分形上颜色是连续变化的；
在分形和非分形的交界处（即 \(A\) 和 \(A^c\) 的边界上）颜色应该是不连续的，从而产生泾渭分明的效果。

做到这一点并不难，但是需要在每一次迭代时考虑当前点的位置信息，这就是所谓的轨道着色。绝大多数惊艳的分形作品都是通过轨道着色技巧来上色的。

我们首先取一个底色，比如说 color=vec3(0)，在每一次迭代中，根据当前位置 p 生成一个颜色，并以一定的权重与 color 混合。理论上颜色的 rgb 的取值范围应该是 \([0,1]\)，但是多数情况下我们要放宽到 \([-1,1]\) 之间，即颜色可以增加也可以减少。否则如果颜色只增不减的话，那么多次迭代以后 rgb 值很可能会溢出，变成白色。此外，随着迭代次数 \(n\) 的增加，\(F^n(K)\) 越来越接近真实的分形 \(A\)，后面加入的颜色的权重应该单调下降，以保证突出分形的细节。这也符合我们的生活直觉：想象一下，当一位画家作画时，在开始的时候他可以浓墨重彩地画一个轮廓，但是越到后面描绘更加精细的部分时，他就会换用更细的画笔，小心地蘸一点颜料。

根据 p 生成颜色的着色方案无穷无尽，请随便发挥你的创造力。一般来说你需要反复试验各种不同的方案才能找到最合适的。下面的例子使用了一种非常流行的染色方案，它以 cos(vec3(0,1,2)) 作为底色，并根据当前时间 iTime 以及坐标 uv 进行调整：

着色器编程基础

我简单介绍一下着色器编程的基本概念。打开 shadertoy 网站，点击右上角的新建按钮，你会看到一个最简单的动画：

左边的窗口是画布，显示渲染的结果；右边窗口是代码编辑器，你在这里书写着色器代码。写完以后，点击编辑器界面左下角的三角形（或者按下 Alt + Enter）查看编译后的效果。

画布是由若干像素组成的，你需要根据每个像素的位置，即它的 fragCoord 值指定一个颜色。这个过程是在 mainImage 函数中实现的：

1	`void mainImage(out vec4 fragColor, in vec2 fragCoord);`

其中 fragCoord 是像素的位置，fragColor 是需要设置的像素颜色。

现实生活中有一个很形象的例子可以帮你理解着色器编程：假设你是一场方阵表演的导演，所有演员排成一个 \(W\times H\) 的方阵，每个演员可以改变自己衣服的颜色。你的任务是发出合适的指令让每个演员根据自己的位置计算出正确的颜色，使得整个方阵呈现出漂亮的图案。

如果你一个一个地对每个演员下指令，张三你应该显示红色，李四你应该显示蓝色，等等 … 对成千上万个演员，这么挨个下指令还不得把人累死？正确的做法是，你应该同时对所有演员发出相同的指令，比如：“每个人，计算自己和方阵中心的距离，小于 10 米的显示红色，大于等于 10 米的显示黑色”。由于每个演员都是一个单独的 GPU 计算单元，他们可以根据你发出的同一条指令，在极短的时间内（毫秒级）分别计算出各自的颜色。计算完毕后，观众应该会看到一个红色的圆。这种基于相同指令并行计算的工作方式就是 GPU 流水线的机制。

我们来试试在着色器中实际编程绘制这个圆的例子。为了方便起见，我们假设方阵的中心是原点，并尝试绘制一个以原点为中心，半径是 0.5 的圆。这只需要三行代码：

首先，每个演员会根据自己在方阵中的位置 fragCoord，计算自己的归一化坐标 p，使得 p 的纵坐标 p.y 位于 \([-1,1]\) 中；
然后，每个演员计算自己到以原点为中心、半径为 0.5 的圆的距离 d。位于圆内部的演员得出的 d 是负数，圆外部的演员得出的 d 是正数；
最终每个人根据 d 的符号来确定自己的颜色。

你可以修改上面代码中的半径，颜色等参数，观察左边画布的变化来体会效果。

这个例子其实蕴含了 shader 编程的一个非常核心的概念，即距离场 (distance field)。当我们想绘制某个图案时，我们通过计算像素到这个图案的距离来对像素进行着色。在绘制 IFS 分形时，我们同样需要计算像素到分形的距离，并根据距离值来确定像素的颜色。这一点在下节会详细介绍。

空间折叠

上面讲到，在着色器里面画 IFS 就是给每个像素指定一个颜色，这个颜色应该由像素对应的点到分形的距离来决定。假设像素对应空间中的点 \(p\)，初始紧集是 \(K\)，我们用迭代 \(n\) 次的结果 \(F^n(K)\) 作为分形的近似，这里 \(n\) 是某个常数，在绝大多数场景下 \(n=30\) 就足够了。于是我们需要计算距离 \(d(p,F^n(K))\) 并根据这个距离值对像素染色。然而，直接计算 \(d(p,F^n(K))\) 是不可行的，因为如果有 \(N\) 个不同的压缩映射 \(f_1,\ldots,f_N\)，那么每次迭代后集合的数目会乘以 \(N\)，这是指数增长的，10 次迭代后就会产生多达 \(10^N\) 个不同的集合。维护如此数量的集合会轻易耗尽内存。这该怎么办呢？

有个巧妙的做法可以绕过这个困难：我们把 \(d(p,F^n(K))\) 中的 \(F^n\) 挪到另一侧并取逆，转而计算 \(d(F^{-n}(p), K)\)！实际上，如果 \(F\) 只包含旋转、平移、反射、缩放、球反演这些变换的话，\(d(p,F^n(K))\) 和 \(d(F^{-n}(p), K)\) 之间存在非常简单的关系，我们可以通过计算后者来得到前者！这个关系的推导我放在注释中介绍。

于是，在着色器编程时，我们需要将压缩迭代映射的步骤倒过来，采取相反的操作：即将逆映射 \(F^{-1}\) 迭代作用在 \(p\) 上，执行足够的迭代次数后，通过计算距离 \(d(F^{-n}(p), K)\) 来给 \(p\) 对应的像素上色。由于 \(F\) 是“一对多” 的映射，所以 \(F^{-1}\) 是 “多对一”的，反复迭代应用 \(F^{-1}\) 会把空间“折叠”。我们实际上是在这个折叠后的空间上作画。

总而言之，在着色器编程中我们真正需要的函数不是 \(F\)，而是 \(F^{-1}\)！

这个先折叠后画图的操作，也可以用一个生活中的例子来形象地解释，即剪纸艺术：

在剪纸过程中，首先把纸张反复折叠，然后在折叠后的纸张上画出某个特定的图形，沿着这个图形裁剪，再将纸张展开得到的就是美丽的图案。展开纸张的操作对应迭代映射 F，它把一个初始的紧集铺开到空间中变成分形；折叠纸张的操作对应 F^{-1}，它把分形折叠回最初的紧集。 — 在剪纸过程中，首先把纸张反复折叠，然后在折叠后的纸张上画出某个特定的图形，沿着这个图形裁剪，再将纸张展开得到的就是美丽的图案。展开纸张的操作对应迭代映射 \(F\)，它把一个初始的紧集铺开到空间中变成分形；折叠纸张的操作对应 \(F^{-1}\)，它把分形折叠回最初的紧集。

在设置 \(F^{-1}\) 的时候，我们完全不必关心每个 \(f_i\) 是什么，我们需要的只是准备一些折叠函数，然后将它们组合起来得到 \(F^{-1}\)。

如果你去看那些分形作品的代码的话，会发现它们几乎都在 \(F^{-1}\) 中使用了 abs 函数。这个函数是最简单的空间折叠函数，它会把整个空间折叠到第一象限。进一步，再叠加关于平面/球面的反射可以产生出更复杂的折叠。

在下面的例子中，我们首先用 p = abs(p) 将整个空间折叠到第一象限，这相当于折纸的时候将纸对折两次；然后只要在第一象限中剪出一个圆，就可以同时在其它象限自动得到另外三个圆：

\(d(p,F^n(K))\) 和 \(d(F^{-n}(p), K)\) 之间的关系：

如果 \(F\) 是平移、旋转、反射这样的保持 Euclidean 距离不变的刚体运动，那么自然有 \[d(p, F(K)) = d(F^{-1}(p), K).\]

如果 \(F\) 是一个缩放变换，比如 \(F(x) = x/s\,(s>1)\)，那么 \[d(p, F(K)) = d(p, 1/s\cdot K) = 1/s\cdot d(s\cdot p, K) = 1/s\cdot d(F^{-1}(p), K).\] 即我们要对折叠以后算出来的距离值再除以 \(s\)。

既然每次迭代 \(F\) 以比例 \(1/s\) 缩小，所以用 \(F^{-1}\) 迭代 \(n\) 次以后累积放大的比例就是 \(s^n\)，我们要将 \(d(F^{-n}(p), K)\) 再除以 \(s^n\) 才是最终正确的距离值。

对球的反演变换，甚至更一般的变换，我们可以通过计算 \(F^{-1}\) 在 \(p\) 处 Jacobian 矩阵的行列式的绝对值，作为 \(p\) 处缩放的近似。把迭代过程中所有这些行列式的绝对值相乘，并对最后得到的 \(d(F^{-n}(p), K)\) 再除以这个乘积，就可以作为 \(d(p, F(K))\) 的近似。

实战演示

我以 Shadertoy 上一个混淆过的作品为例子来完整展示上面的理论。下面是重新改写后的可读版本，我在注释中解释了每一步的含义：

下一篇：咖啡杯中的焦散线