## 使用 WebGL 提升可视化中的布局性能

文 / 沈毅

现在 PC 端和移动端的浏览器对于 WebGL 的支持都已非常普及，我们常常会利用 WebGL 强大的特性去绘制三维的场景用于可视化、广告展示、产品模型展示和小的休闲游戏。除了三维渲染这个 WebGL 最擅长的老本行，大家也会尝试使用 WebGL 去替代 Canvas 或者 SVG，加速二维图形的渲染，GIS 服务商 Mapbox 就在提供的 JS SDK 中使用 WebGL 来绘制地图。

除了三维或者二维的渲染，其实我们也可以利用 WebGL 来做 GPU 通用计算，去加速一些适合并行的算法，例如物理计算、布料模拟、FFT 等等。这篇文章介绍我们是如何在 ECharts GL 中使用 WebGL 对力引导布局提升近百倍的性能的。

### 什么是力引导算法

使用过 ECharts 的同学可能知道 ECharts 中有个系列类型是 graph，这个系列类型可以用来可视化表示节点—边的关系数据或者网络数据。而用户在刚拿到关系数据或者网络数据的时候是没有节点位置等布局信息的，大部分都只有每个节点的权重、类别等信息，以及表示两个关系的边。如何从这些有限的信息布局得到美观的网络图一直是网络关系数据的可视化中的研究重点。

关系图中最常用的布局算法是力引导（Force Directed Layout）算法，力引导算法的思想非常朴素，算法将关系数据中的每个节点模拟为电荷，将关系数据中的每条边模拟为弹簧，这样任意两个节点之间存在一个电荷上的斥力，如果两个节点之间存在关系边，那这两个节点还会受到一个弹簧的引力，所以力引导算法也会叫做 Spring Embedded Layout。

每次迭代的时候，我们依次计算每个节点来自其它节点的斥力以及节点的邻接边所产生的引力的总和，然后根据这个受力对节点位移。第二次迭代的时候就根据新的位置继续计算受力和位移，就这样一直迭代到所有节点都受力平衡了，整个布局的过程就算完成了。这个时候整个图的能量处于最小值，节点之间关系边因为互相的引力会聚合在一起，而不同的圈子因为节点之间的斥力会分散开来，因此可以形成不同圈子聚类的形态。

图1就是对一份典型的微博转发关系的网络数据进行力引导布局的结果。

<img src="http://ipad-cms.csdn.net/cms/attachment/201707/5954b63a90faf.png" alt="图1  微博转发关系图" title="图1  微博转发关系图" />

图1  微博转发关系图

### 力引导算法的性能问题和常见优化

可以看到力引导算法中的每次迭代都需要两两计算节点之间的受力，所以每次迭代的开销都是O(n2)，而且往往整个布局需要几百次迭代后才能呈现出比较好的形态，导致整个布局的性能开销往往非常大。为了防止长时间布局带来的 UI 阻塞，很多时候我们会将布局的过程动态地呈现出来，也就是每次节点位置发生变化都会重绘更新画面。并且通过 JavaScript 中的 requestAnimationFrame 接口将布局和绘制的迭代分摊到不同帧里，这样既有效地防止了布局对 UI 的阻塞，同时因为算法是基于物理模拟的，所以整个布局的动画往往非常有意思，也增加了这个布局算法的趣味性。

除了单次迭代的开销比较大之外，传统的力引导算法对稍微大点规模的关系图也存在着布局收敛速度过慢的问题，因为节点的位移是离散的，所以很多节点会在受力均衡点摆动（Swing）而无法到达受力均衡点，会导致布局长时间无法结束。

大部分力引导算法的论文对大规模网络图布局在性能上的改进都集中在上面两点上，减小单次迭代的开销以及加速布局收敛的速度。

第一点减小单次迭代的开销，最常见的就是使用 Barnes Hut Simulation 计算节点之间的受力。这个方法将整个空间用四叉树划分，每个关系图的节点会放入四叉树的叶子节点中，这样距离比较远的两片区域可以不用每个节点之间都计算斥力了，只需要将两片区域作为整体去计算斥力。理论上使用 Barnes Hut Simulation 后可以从原来 O(n2) 的时间复杂度降低到 O(n log n) 的时间复杂度。

第二点加快整个布局的收敛速度，也就是减少布局的迭代次数。最早在《Graph Drawing by Force-directed Placement》一文中提出了使用模拟退火（Simulated Annealing）加速整个布局的稳定和收敛，模拟退火引入了一个温度（Temporature）因子，跟物理意义一样，温度会影响每个节点的活跃度，也就是影响节点受力后位移的距离，布局刚开始的时候温度比较高，节点受力后移动的范围比较大，能够快速接近受力平衡的位置，随着布局的进行，温度会逐渐降低，节点的移动幅度也越来越小，保证节点不会因为单次移动幅度过大而产生在平衡位置附近摇摆的情况。

除了算法层面的优化，在程序层面，我们可以利用 WebWorker 将布局放到一个单独的线程里，跟绘制渲染的主线程分开来，布局不会阻塞绘制，绘制也不会阻塞掉下一次的布局。既可以加快整体布局的速度，也可以让整个布局的动画更流畅。

这些优化能够让一个上千节点的关系图从原来的上百 ms 一次迭代优化到10ms一次迭代，布局的动画也更加流畅。但是我们依然没法做到短时间内就让整个布局完成。一方面是因为算法本身复杂度的限制，还有一方面是 JavaScript 多线程的限制导致无法发挥 CPU 多核的性能，也无法利用 SIMD 去加速向量的计算。 

在算法复杂度上现在基本上已经到达瓶颈了，那么在程序的实现层面，我们是不是可以利用 WebGL 实现力引导布局的 GPU 通用计算（GPGPU）从而绕过 JavaScript 的限制。而且因为力引导算法每个节点的受力计算互相不受影响，所以节点的受力很适合在 GPU 中并行计算。

### WebGL 中实现 GPU 通用计算

很早之前就有人开始尝试在 WebGL 中通过 GPGPU 去加速布料模拟（Cloth Simulation），群体模拟（Flocking），粒子系统（Particle System）等适合 GPU 加速的算法，并且取得了还不错的效果。例如 three.js 中就有分别用 CPU 和 GPU 实现 Flocking 的例子，两者的帧率在群体数量为4096的时候可以相差几十倍，如图2所示。

<img src="http://ipad-cms.csdn.net/cms/attachment/201707/5954b724c7f9e.png" alt="图2  帧率上的较大差异" title="图2  帧率上的较大差异" />

图2  帧率上的较大差异

shader 是 WebGL 里的可编程着色器，我们通常是用 shader 来计算顶点的位置，像素的光照等三维渲染方面的东西。如果要直接使用 WebGL 1.0 的 shader 做 GPU 的通用计算，有个很大的限制是 shader 无法直接读写显存，shader 的输入只有纹理数据和顶点属性（Attributes）的数据。而 shader 最终的输出也只有写入到 Drawing Buffer 的 RGBA 颜色。

所以通常在 WebGL 中实现通用计算的思路是将需要计算的数据存储在纹理里，然后在着色器（Shader）里读取浮点纹理中的数据，做一定的运算后再配合 Framebuffer 写入到另一个纹理中，然后将两个纹理交换一下，现在写入的纹理作为下一次读的纹理，也就是说将纹理作为暂存数据的地方，见图3所示。随着 WebGL 中浮点纹理扩展`GL_ARB_texture_float`的普及，RGBA 每个通道能存储的数据精度也从原来的8字节变成了32字节，对于大部分数据来说这个精度已经够用了。

<img src="http://ipad-cms.csdn.net/cms/attachment/201707/5954b749e269d.png" alt="图3  GPU通用计算的流程" title="图3  GPU通用计算的流程" />

图3  GPU 通用计算的流程

下面是根据输入的存储有节点受力和之前位置的纹理，计算得到新位置的 shader 代码。

```
uniform sampler2D forceTex;
uniform sampler2D forcePrevTex;
uniform sampler2D positionTex;
uniform sampler2D globalSpeedTex;
varying vec2 v_Texcoord;
void main() {
    // 从纹理中获取节点的受力
    vec2 force = texture2D(forceTex, v_Texcoord).xy;
    vec2 forcePrev = texture2D(forcePrevTex, v_Texcoord).xy;
    // 从纹理中获取节点之前的位置
    vec4 node = texture2D(positionTex, v_Texcoord);
 
    // 计算节点需要移动的距离
    float globalSpeed = texture2D(globalSpeedTex, vec2(0.5)).r;
    float swing = length(force - forcePrev);
    float speed = 0.1 * globalSpeed / (0.1 + globalSpeed * sqrt(swing));
 
    float df = length(force);
    if (df > 0.0) {
        speed = min(df * speed, 10.0) / df;
        // 根据之前的位置，需要移动的距离，计算得到新的位置后输出
        gl_FragColor = vec4(node.xy + speed * force, node.zw);
    }
    else {
        gl_FragColor = node;
    }
}
```

### 在 WebGL 中实现力引导布局

我们可以从刚才那段代码中一窥怎么在 shader 中做位置计算，实际上这个位置的更新是力引导布局中每次迭代计算的最后一步，也是实现相对简单的一步，我们接下来就来详细讲讲整个力引导布局如何在 WebGL 中实现。

#### 初始化节点和边的数据

力引导算法会先基于某个简单高效的算法初始化每个节点的位置，合理的初始化方式可以加速后面布局的收敛，但是一般情况下就是在一个范围内随机分布节点。

节点位置初始化完之后我们需要将其存入到浮点纹理中作为初始的来源。WebGL 除了支持加载的图片作为纹理数据之外，也支持使用 TypedArray 作为纹理像素数据。浮点纹理的话则是 Float32Array。

```
var offset = 0; 
for (var i = 0; i < nodes.length; i++) {
    var node = nodes[i];
    positionBuffer[offset++] = node.x;
    positionBuffer[offset++] = node.y;
    positionBuffer[offset++] = node.mass;
    positionBuffer[offset++] = node.size;
}
```

一张纹理总共有四个通道 RGBA 可以利用，我们在 RG 通道中存储节点的 x，y 位置，在 B 通道中存储节点的质量/权重，在 A 通道中存储节点的半径，这个半径计算节点之间是否重叠的时候可以派上用场。

除了节点的数据，我们还需要初始化边的数据。对于节点来说，能够做到输入的节点像素和输出像素一一对应，但是对于边来说不是，我们在计算边对两个节点的引力的时候，需要能够在 shader 里正确索引到相应的节点，并且计算完之后能够写入正确的节点。

因此跟节点信息存在纹理里不一样，我们把边的两个节点和权重存在了 WebGL 的 attributes 中。

```
for (var i = 0; i < edges.length; i++) {
    var attributes = edgeGeometry.attributes;
    var weight = edges[i].weight;
    if (weight == null) {
        weight = 1;
    }
    attributes.node1.set(i, this.getNodeUV(edges[i].node1, uv));
    attributes.node2.set(i, this.getNodeUV(edges[i].node2, uv));
    attributes.weight.set(i, weight);
}
```

其中 getNodeUV 方法将节点的线性索引映射到了二维的纹理坐标，方便在 shader 中根据这个纹理坐标索引到正确的节点。

这些放入显存中的数据在接下来每次迭代计算受力的时候会被频繁的用到。

#### 计算受力

初始化完之后就是每次迭代计算节点与节点的斥力以及边的引力。

图4就是每次迭代的流程，整个流程分为三步：1. 计算节点与节点之间的斥力后写入存储每个节点受力的 Force 纹理；2. 计算边对两个节点的引力后叠加到 Force 纹理；3. 根据 Force 纹理计算得到每个节点新的位置写入 Position 纹理。

<img src="http://ipad-cms.csdn.net/cms/attachment/201707/5954bdbe736aa.png" alt="图4  每次迭代的流程" title="图4  每次迭代的流程" />

图4  每次迭代的流程

在计算节点与节点之间的斥力时，我们用了一个全屏的 Pass（Fullscreen Quad Pass）渲染到一个 FrameBuffer 所绑定的纹理对象，这个纹理对象的每个像素对应的就是每个节点所受到的力。所以这一步是在像素着色器（Fragment Shader）里做的。

```
#define NODE_COUNT 0
// 存储节点位置等信息的纹理 
uniform sampler2D positionTex;
 
uniform vec2 textureSize;
uniform float scaling;
uniform float gravity;
uniform vec2 gravityCenter;
 
varying vec2 v_Texcoord;
 
void main() {
 
    vec4 n0 = texture2D(positionTex, v_Texcoord);
 
vec2 force = vec2(0.0);
    // 遍历所有其它节点计算受力总和
for (int i = 0; i < NODE_COUNT; i++) {
        // 根据节点的索引计算节点的纹理坐标
        vec2 uv = vec2(
            mod(float(i), textureSize.x) / (textureSize.x - 1.0),
            floor(float(i) / textureSize.x) / (textureSize.y - 1.0)
        );
        // 从纹理中获取节点信息
        vec4 n1 = texture2D(positionTex, uv);
 
        vec2 dir = n0.xy - n1.xy;
        float d2 = dot(dir, dir);
 
        if (d2 > 0.0) {
            float factor = scaling * n0.z * n1.z / d2;
            force += dir * factor;
        }
}
 
    // 计算向心力
    vec2 dir = gravityCenter - n0.xy;
    force += dir * n0.z * gravity;
 
    gl_FragColor = vec4(force, 0.0, 1.0); 
}
```

一些类似向心力的因子，布局的缩放因子等全局的布局参数会通过 uniform 传入给 shader 使用。

接下来是计算边的引力。因为要绘制到同一个存储有节点受力的 Force 纹理中而且要计算受力的叠加，我们需要开启 WebGL 的混合模式并把混合模式设置成叠加。

```
gl.enable(gl.BLEND);
gl.blendEquation(gl.FUNC_ADD);
gl.blendFunc(gl.ONE, gl.ONE);
```

然后在顶点着色器中去计算边的引力。

```
// 边的第一个节点的纹理坐标
attribute vec2 node1;
// 边的第二个节点的纹理坐标
attribute vec2 node2;
attribute float weight;
 
// 存储节点位置等信息的纹理
uniform sampler2D positionTex;
     
uniform float edgeWeightInfluence;
uniform vec2 windowSize;
 
varying vec2 v_Force;
 
void main() {
 
    vec4 n0 = texture2D(positionTex, node1);
    vec4 n1 = texture2D(positionTex, node2);
 
    vec2 dir = n1.xy - n0.xy;
    float d = length(dir);
float w = pow(weight, edgeWeightInfluence);
    // 计算这次计算需要写入的像素位置，也就是存储第一个节点受力的像素。
    vec2 offset = vec2(1.0 / windowSize.x, 1.0 / windowSize.y);
    vec2 scale = vec2((windowSize.x - 1.0) / windowSize.x, (windowSize.y - 1.0) / windowSize.y);
    vec2 pos = node1 * scale * 2.0 - 1.0;
gl_Position = vec4(pos + offset, 0.0, 1.0);
// 使用 gl.POINT 的画点模式
    gl_PointSize = 1.0;
    if (d <= 0.0) {
        v_Force = vec2(0.0);
        return;
}
// v_Force 会传给像素着色器后直接写入到 gl_FragColor 中
    v_Force = dir * w;
}
```

计算完受力后，就可以根据计算得到这张存储中每个节点受到的总的力的纹理，进一步计算出节点位移的速度以及根据这个速度位移得到的新的位置了。

#### 判断布局结束

我们判断布局是否结束主要是通过计算所有节点总的受力，如果受力小于某个阈值则可以判断整个布局是稳定下来了。但是可以从前面的实现中看到，每一次迭代，JavaScript 除了负责整个布局迭代的调度，调用 WebGL 的绘制接口，设置绘制参数，基本上不做任何布局计算的事。所有的布局计算都是在 GPU 中完成，布局的中间数据和结果数据，包括节点受力也是存在显存中。所以每一次迭代都需要 JavaScript 从显存中回读节点的数据，然后判断布局是否结束，结束的话就不通过 requestAnimationFrame 启用下一次迭代了。

WebGL 提供了 readPixels 方法从显存的 Drawing Buffers 中读取绘制出来的像素数据。

```
var forceArr = new Float32Array(width * height * 4);
gl.readPixels(
    0, 0, width, height, gl.RGBA, gl.FLOAT, forceArr
);
```

需要注意的是，目前只有 Chrome 浏览器支持对浮点纹理的回读，其它浏览器只支持读取`gl.UNSIGNED_BYTE`格式的像素数据。所以如果要支持其它浏览器，例如 Safari，我们需要多一个绘制的 Pass 将所有节点的受力总和存储到一张1x1大小，只有一个像素的`UNSIGNED_BYTE`格式的纹理中。这个求和过程跟计算边的受力是一样的，就是在开启叠加混合的模式后，将所有节点的力都输出到一个像素中。

```
var forceSum = new Float32Array(4);
gl.readPixels(0, 0, 1, 1, gl.RGBA, gl.UNSIGNED_BYTE, forceSum);
```

因为我们只需要一个是否超过阈值的判断，所以`UNSIGNED_BYTE`这个精度对于我们来说够用了。

### 性能对比

那么利用 WebGL 加速后性能上究竟能有多少提升。作为测试，我们挑了一份 Gephi 提供的比较大的关系图数据，见图5。这份关系数据拥有大约22k的节点和48k的边。分别在一台2012年的13寸 Macbook Pro 和一台最近的 GTX 1070显卡+ i7 CPU 的台式机上做了测试。

<img src="http://ipad-cms.csdn.net/cms/attachment/201707/5954bfb1ced9d.png" alt="图5  拥有22k节点，48k边的关系图" title="图5  拥有22k节点，48k边的关系图" />

图5  拥有22k节点，48k边的关系图

#### 在13寸 Macbook Pro 上一次布局迭代：

JavaScript：

- without Barnes Hut：~28000ms
- with Barnes Hut：~1000ms

WebGL：~260ms

#### 在 GTX 1070 + i7的台式机上一次布局迭代：

JavaScript：

- without Barnes Hut：~12000ms
- with Barnes Hut：~300ms

GPU：~2ms

可以看到，相对于 JavaScript 的实现，WebGL 的实现对于高端显卡的提速非常明显，甚至已经到达了上百倍的性能提升。

### 总结

这篇文章我们介绍了 ECharts 中对于关系图的力引导布局的算法和常见的性能优化方案，以及使用 WebGL 实现力引导布局的思路，包括如何计算节点的引力和斥力，如何判断布局是否结束等。

WebGL 对于力引导布局这种适合并行的算法有着非常可观的性能提升，特别是电脑使用的是高端显卡的时候，一方面是因为现代显卡的显存带宽和浮点运算非常可观，另一方面因为 JavaScript 本身的限制导致只能利用 CPU 的部分性能。

而且现在支持 WebGL 和浮点纹理的浏览器也已经十分普及了，我们除了能够用 WebGL 来绘制一些三维的场景和酷炫的特效之外，完全可以尝试一下利用 WebGL 去做一些适合并行加速的通用计算，说不定能收获到出乎意料的性能提升。